1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
import numpy as np
import pandas as pd
from statsmodels.tsa.stattools import grangercausalitytests

# 예시 데이터 생성
np.random.seed(1)
x = np.random.randn(1000)
y = np.random.randn(1000)
y[1:] += 0.5 * x[:-1]  # X가 Y를 그랜저 인과하도록 설정

data = pd.DataFrame({'X': x, 'Y': y})

# 그랜저 인과관계 검정
result = grangercausalitytests(data[['Y', 'X']], maxlag=5)

# 결과 출력
for lag, values in result.items():
    print(f"Lag {lag}:")
    print(f"  F-statistic: {values[0]['ssr_ftest'][0]}")
    print(f"  p-value: {values[0]['ssr_ftest'][1]}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
from statsmodels.tsa.api import VAR
from statsmodels.tsa.stattools import adfuller

# 데이터 정상화 (필요한 경우)
def difference(data):
    return data.diff().dropna()

# ADF 테스트로 정상성 확인
def adf_test(series):
    result = adfuller(series)
    print(f'ADF Statistic: {result[0]}')
    print(f'p-value: {result[1]}')

# VAR 모델 적용
model = VAR(data)
results = model.fit(maxlags=15, ic='aic')

# 그랜저 인과관계 검정
granger_results = results.test_causality('Y', 'X', kind='f')
print(granger_results.summary())

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
from statsmodels.tsa.arima.model import ARIMA
import matplotlib.pyplot as plt

# ARIMA 모델 적용
model = ARIMA(data['Y'], order=(1,1,1))
results = model.fit()

# 예측
forecast = results.forecast(steps=10)

# 결과 시각화
plt.figure(figsize=(12,6))
plt.plot(data['Y'], label='Observed')
plt.plot(forecast, label='Forecast')
plt.legend()
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 필요 라이브러리 로드
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

import warnings
warnings.filterwarnings('ignore')

from statsmodels.tsa.api import VAR
import statsmodels.tsa.api as smt
from statsmodels.tsa.stattools import grangercausalitytests
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.arima.model import ARIMA

1
2
3
4
5
6
# 데이터 로드 (데이터는 직접 다운 받아주세요)
# https://www.kaggle.com/datasets/yasserh/advertising-sales-dataset
df = pd.read_csv("./Data/Advertising Budget and Sales.csv")

print(df.shape)
df.head()

1
2
3
4
5
6
7
variables = ['TV Ad Budget ($)', 'Radio Ad Budget ($)', 'Newspaper Ad Budget ($)', 'Sales ($)']

# 각 값의 첫 값으로 정규화
df['Sales ($)'] = df['Sales ($)'] / df['Sales ($)'].iloc[0]
df['TV Ad Budget ($)'] = df['TV Ad Budget ($)'] / df['TV Ad Budget ($)'].iloc[0]
df['Radio Ad Budget ($)'] = df['Radio Ad Budget ($)'] / df['Radio Ad Budget ($)'].iloc[0]
df['Newspaper Ad Budget ($)'] = df['Newspaper Ad Budget ($)'] / df['Newspaper Ad Budget ($)'].iloc[0]

1
2
df.plot(figsize=(15,3))
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# matplotlib의 기본 색상 순서 사용
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

# 각 변수를 별도의 서브플롯으로 그리기
fig, axes = plt.subplots(len(variables), 1, figsize=(15, 3*len(variables)))
fig.suptitle('Normalized Values Over Time (Subplots)', fontsize=16)

for i, var in enumerate(variables):
    axes[i].plot(df.index, df[var], label=var, color=colors[i])
    axes[i].set_title(var)
    axes[i].set_xlabel('Index')
    axes[i].set_ylabel('Normalized Value')
    axes[i].legend()
    axes[i].grid(True)

    # y축 범위 설정 (0부터 데이터의 최대값까지)
    axes[i].set_ylim(0, df[var].max() * 1.1)  # 최대값에 10% 여유 추가

plt.tight_layout()
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
from statsmodels.tsa.stattools import adfuller, grangercausalitytests

# ADF 테스트 함수 정의
def adf_test(timeseries):
    result = adfuller(timeseries, autolag='AIC')
    return pd.Series({'Test Statistic': result[0], 'p-value': result[1], 'Critical Values': result[4]})

adf_results = {var: adf_test(df[var]) for var in variables}

# ADF 테스트 결과 출력 및 해석
for var, result in adf_results.items():
    print(f"\nADF Test Results for {var}:")
    print(result)
    
    if result['p-value'] <= 0.05:
        print(f"해석: {var}는 5% 유의수준에서 정상 시계열입니다.")
    else:
        print(f"해석: {var}는 5% 유의수준에서 비정상 시계열입니다.")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
# 결과
ADF Test Results for TV Ad Budget ($):
Test Statistic                                            -14.158141
p-value                                                          0.0
Critical Values    {'1%': -3.4636447617687436, '5%': -2.876176117...
dtype: object
해석: TV Ad Budget ($)는 5% 유의수준에서 정상 시계열입니다.

ADF Test Results for Radio Ad Budget ($):
Test Statistic                                            -14.129897
p-value                                                          0.0
Critical Values    {'1%': -3.4636447617687436, '5%': -2.876176117...
dtype: object
해석: Radio Ad Budget ($)는 5% 유의수준에서 정상 시계열입니다.

ADF Test Results for Newspaper Ad Budget ($):
Test Statistic                                            -13.344195
p-value                                                          0.0
Critical Values    {'1%': -3.4636447617687436, '5%': -2.876176117...
dtype: object
해석: Newspaper Ad Budget ($)는 5% 유의수준에서 정상 시계열입니다.

ADF Test Results for Sales ($):
Test Statistic                                            -13.990162
p-value                                                          0.0
Critical Values    {'1%': -3.4636447617687436, '5%': -2.876176117...
dtype: object
해석: Sales ($)는 5% 유의수준에서 정상 시계열입니다.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import grangercausalitytests
from statsmodels.tsa.api import VAR
from statsmodels.tools.eval_measures import aic, bic

def optimal_lag(data, max_lag):
    model = VAR(data)
    results = {}
    for lag in range(1, max_lag + 1):
        result = model.fit(lag)
        results[lag] = result.aic
    return min(results, key=results.get)

def granger_causality(data, variables, max_lag=12, alpha=0.05):
    results = {}
    for i in range(len(variables)):
        for j in range(len(variables)):
            if i != j:
                pair = (variables[i], variables[j])
                opt_lag = optimal_lag(data[[variables[i], variables[j]]], max_lag)
                granger_result = grangercausalitytests(data[[variables[i], variables[j]]], maxlag=opt_lag, verbose=False)
                p_values = [granger_result[lag][0]['ssr_ftest'][1] for lag in range(1, opt_lag+1)]
                
                results[pair] = {
                    'optimal_lag': opt_lag,
                    'p_values': p_values,
                    'significant': any(p < alpha for p in p_values)
                }
                
                plt.figure(figsize=(10, 6))
                plt.plot(range(1, opt_lag+1), p_values, marker='o')
                plt.axhline(y=alpha, color='r', linestyle='--')
                plt.title(f'Granger Causality: {variables[i]} -> {variables[j]}')
                plt.xlabel('Lag')
                plt.ylabel('p-value')
                plt.show()
    
    return results

# 그랜저 인과관계 분석 수행
results = granger_causality(df, variables)

# 결과 출력
for pair, result in results.items():
    print(f"\nGranger Causality: {pair[0]} -> {pair[1]}")
    print(f"Optimal lag: {result['optimal_lag']}")
    print(f"Significant: {result['significant']}")
    print(f"p-values: {result['p_values']}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.api import VAR
from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_squared_error
from math import sqrt

# 데이터 로드
df = pd.read_csv("./Data/Advertising Budget and Sales.csv")
df = df[['TV Ad Budget ($)', 'Sales ($)']]

df['Sales ($)'] = df['Sales ($)'] / df['Sales ($)'].iloc[0]
df['TV Ad Budget ($)'] = df['TV Ad Budget ($)'] / df['TV Ad Budget ($)'].iloc[0]

df.columns = ['tv', 'sales'] # 간단하게 이름 변경

# 데이터 분할 (80% 훈련, 20% 테스트)
train_size = int(len(df) * 0.8)
train, test = df[:train_size], df[train_size:]

# 그래프 설정
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10), sharex=True)
fig.suptitle('Train and Test Data Visualization')

# TV Ad Budget 그래프
ax1.plot(df.index[:train_size], train['tv'], label='Train')
ax1.plot(df.index[train_size:], test['tv'], label='Test')
ax1.axvline(x=train_size, color='r', linestyle='--', label='Train/Test Split')
ax1.set_ylabel('TV Ad Budget ($)')
ax1.legend()

# Sales 그래프
ax2.plot(df.index[:train_size], train['sales'], label='Train')
ax2.plot(df.index[train_size:], test['sales'], label='Test')
ax2.axvline(x=train_size, color='r', linestyle='--', label='Train/Test Split')
ax2.set_ylabel('Sales ($)')
ax2.legend()

# x축 레이블 설정
ax2.set_xlabel('Time')

plt.tight_layout()
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# USING VAR
lag = 2

# VAR 모델
model_var = VAR(train)
results_var = model_var.fit(lag)

# VAR 동적 예측
forecast_var = []
forecast_input_var = train.values[-lag:]

for i in range(len(test)):
    fc = results_var.forecast(forecast_input_var, steps=1)
    forecast_var.append(fc[0])
    
    # 실제 test 데이터의 이전 lag 값으로 forecast_input_var 업데이트
    if i + 1 < len(test):
        forecast_input_var = np.vstack([forecast_input_var[1:], test.iloc[i].values])
    
forecast_var = np.array(forecast_var)
rmse_var_sales = sqrt(mean_squared_error(test['sales'], forecast_var[:, 1]))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
train_sales = train['sales']
test_sales = test['sales']

# ARIMA 모델
model_arima = ARIMA(train_sales, order=(lag, 0, 0))
results_arima = model_arima.fit()

# 동적 예측을 위한 ARIMA
forecast_arima_dynamic = []
history = train_sales.tolist()

for t in range(len(test_sales)):
    model_arima_dynamic = ARIMA(history, order=(lag, 0, 0))
    model_fit_dynamic = model_arima_dynamic.fit()
    yhat = model_fit_dynamic.forecast(steps=1)[0]
    forecast_arima_dynamic.append(yhat)
    history.append(test_sales.iloc[t])

rmse_arima_sales_dynamic = sqrt(mean_squared_error(test_sales, forecast_arima_dynamic))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
print("\nModel Comparison:")
print(f"VAR Model RMSE: {rmse_var_sales:.5f}")
print(f"VAR Model MSE: {mse_var_sales:.5f}")
print(f"VAR Model MAE: {mae_var_sales:.5f}")
print("-"*50)
print(f"ARIMA Model RMSE: {rmse_arima_sales_dynamic:.5f}")
print(f"ARIMA Model MSE: {mse_arima_sales_dynamic:.5f}")
print(f"ARIMA Model MAE: {mae_arima_sales_dynamic:.5f}")

# Model Comparison:
# VAR Model RMSE: 0.23431
# VAR Model MSE: 0.05490
# VAR Model MAE: 0.19481
# --------------------------------------------------
# ARIMA Model RMSE: 0.23493
# ARIMA Model MSE: 0.05519
# ARIMA Model MAE: 0.18524

1
2
3
4
5
6
7
8
9
10
11
12
# 예측 결과 시각화
plt.figure(figsize=(12, 6))
plt.plot(test.index, test.sales, label='Actual Sales', color='black')
plt.plot(test.index, test.tv, label='Actual TV', color='orange', linestyle = '--')

plt.plot(test.index, forecast_var[:, 1], label='VAR Sales Forecast', color='blue')
plt.plot(test.index, forecast_arima_dynamic, label='Dynamic ARIMA Sales Forecast', color='red')
plt.legend()
plt.xlabel('Time')
plt.ylabel('Sales ($)')
plt.title('Sales Forecast Comparison: VAR vs Dynamic ARIMA')
plt.show()