import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

data = np.genfromtxt("./data_science/07. web_traffic.tsv")
data

# 첫번째 열 [ : , 0 ]은 시간(x)
# 두번째 열 [ : , 1 ]은 hit수(y) 
x = data[:, 0]
y = data[:, 1]

# 의미없는 데이터 확인
np.sum( np.isnan(y) )

# 사용가능한 데이터 얻기
x = x[~np.isnan(y)]
y = y[~np.isnan(y)]

# 그래프

plt.figure( figsize=(10, 6) )
plt.scatter( x, y, s=10 )

plt.title("Web traffic over the last month")
plt.xlabel("Time")
plt.ylabel("Hits/hour")
plt.grid( True, linestyle='-', color='0.75' )
# 데이터가 가급적 화면에 꽉 차도록 autoscale을 적용함 
plt.autoscale( tight=True )
# plt.xticks([A], [B]) ~ A라고 표기될 x축 라벨을 B로 표기
plt.xticks( [w*7*24  for  w  in  range(5)] , [ 'week %i' % w for w in range(5) ] )
plt.show()

def error(f, x, y):
    return np.sqrt( np.mean( f(x) - y )**2 )

# 1차 함수에서 기울기와 y절편을 구한다. 

np.polyfit( x, y, 1)

# 위 기울기와 y절편으로 1차 함수를 만들기

fp1 = np.polyfit( x, y, 1 )
f1 = np.poly1d( fp1 )

error(f1, x, y)

fx = np.linspace( 0, x[-1], 1000 )

plt.figure( figsize=(10, 6) )
plt.scatter( x, y, s=10 )
# 1차 추쇄선 삽입
plt.plot( fx, f1(fx), lw=4 )

plt.title("Web traffic over the last month")
plt.xlabel("Time")
plt.ylabel("Hits/hour")
plt.grid( True, linestyle='-', color='0.75' )
# 데이터가 가급적 화면에 꽉 차도록 autoscale을 적용함 
plt.autoscale( tight=True )
# plt.xticks([A], [B]) ~ A라고 표기될 x축 라벨을 B로 표기
plt.xticks( [w*7*24  for  w  in  range(5)] , [ 'week %i' % w for w in range(5) ] )
plt.legend( [ "d=%i" % f1.order], loc=2 )
plt.show()

f2p = np.polyfit(x, y, 2)
f2 = np.poly1d( f2p )

f3p = np.polyfit( x, y, 3 )
f3 = np.poly1d( f3p )

f15p = np.polyfit( x, y, 15 )
f15 = np.poly1d(f15p)

error( f1, x, y), error( f2, x, y), error(f3, x, y), error(f15, x, y)

plt.figure( figsize=(10, 6) )
plt.scatter( x, y, s=10 )

# 추쇄선 삽입
plt.plot( fx, f1(fx), lw=4, label='f1')
plt.plot( fx, f2(fx), lw=4, label='f2')
plt.plot( fx, f3(fx), lw=4, label='f3')
plt.plot( fx, f15(fx), lw=4, label='f15')

plt.title("Web traffic over the last month")
plt.xlabel("Time")
plt.ylabel("Hits/hour")
plt.grid( True, linestyle='-', color='0.75' )
# 데이터가 가급적 화면에 꽉 차도록 autoscale을 적용함 
plt.autoscale( tight=True )
# plt.xticks([A], [B]) ~ A라고 표기될 x축 라벨을 B로 표기
plt.xticks( [w*7*24  for  w  in  range(5)] , [ 'week %i' % w for w in range(5) ] )
plt.legend( loc=2 )
plt.show()

import warnings
warnings.filterwarnings("ignore")

import itertools
import pandas as pd
import numpy as np

import statsmodels.api as sm
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
%matplotlib inline

df = sm.datasets.get_rdataset("AirPassengers").data

df.head(3)

def yearfraction2datetime( yearfraction, startyear=0 ):
    import datetime, dateutil
    
    year = int( yearfraction ) + startyear
    month = int( round( 12*( yearfraction - year) ) )
    delta = dateutil.relativedelta.relativedelta( months=month )
    date = datetime.datetime( year, 1, 1) + delta
    return date

yearfraction2datetime(1949.333333)

import datetime, dateutil

dateutil.relativedelta.relativedelta( months=15 )

df["datetime"] = df.time.map( yearfraction2datetime )

df.tail(3)

t = df.time
y = df.AirPassengers

plt.figure(figsize=(12, 8))
plt.plot(t, y)
plt.show()

result = sm.OLS.from_formula( "AirPassengers ~ time", data=df ).fit()
print( result.params )
print( result.params[0], result.params[1] )

t = df.time
y = df.AirPassengers
trend = result.params[0] + result.params[1] * t

plt.figure(figsize=(12, 8))
plt.plot(t, y, 'o-', t, trend, '-')
plt.show()

result2 = sm.OLS.from_formula( "AirPassengers ~ time + np.power(time, 2)", data=df ).fit()
print( result2.params[0], result2.params[1], result2.params[2] )

trend2 = result2.params[0] + result2.params[1] * t + result2.params[2]* t**2
plt.figure(figsize=(12, 8))
plt.plot(t, y, 'o-', t, trend2, '-')
plt.show()

%matplotlib inline
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import datetime
from dateutil.relativedelta import relativedelta
import seaborn as sns
import statsmodels.api as sm
from statsmodels.tsa.stattools import acf
from statsmodels.tsa.stattools import pacf
from statsmodels.tsa.seasonal import seasonal_decompose

df = pd.read_csv( './data_science/07. portland-oregon-average-monthly-.csv', index_col=0 )
df.index.name=None
df.reset_index(inplace=True)
df.head(3)

start = datetime.datetime.strptime("1973-01-01", "%Y-%m-%d")
date_list = [start + relativedelta(months=x) for x in range(0, len(df)) ]
df['index'] = date_list
df.set_index( ['index'], inplace=True) 
df.index.name=None
df.head(3)

df.columns = ['riders']
df['riders'] = df.riders.apply( lambda x: int(x)*100 ) 
df.head()

df.riders.plot(
    figsize=(12, 8),  
    title='Monthly Ridership',
    fontsize=14
)

plt.savefig(
    'month_ridership.png',
    bbox_inches='tight'
)

decomposition = seasonal_decompose( df.riders, freq=12 )
fig = plt.figure()
fig = decomposition.plot()
fig.set_size_inches(12, 8)

trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid

# Define the p, d and q parameters to take any value between 0 and 2
p = d = q = range(0, 2)

# Generate all different combinations of p, q and d triplets
pdq = list( itertools.product(p, d, q))

# Generate all different combinations of seasonal p, q and d triplets
seasonal_pdq = [ (x[0], x[1], x[2], 12)  for x in pdq ]

print('Example of parameter combinations for Seasonal ARIMA ...')
print('SARIMAX:  {} x {}'.format(pdq[1], seasonal_pdq[1]) )
print('SARIMAX:  {} x {}'.format(pdq[1], seasonal_pdq[2]) )
print('SARIMAX:  {} x {}'.format(pdq[2], seasonal_pdq[3]) )
print('SARIMAX:  {} x {}'.format(pdq[2], seasonal_pdq[4]) )

warnings.filterwarnings("ignore")

select_candi = 10000000
param_candi = ( 0, 0, 0 )
param_seasonal_candi = ( 0, 0, 0)

count=0
end_count = len(pdq)

for param in pdq:   
    for param_seasonal in seasonal_pdq:
        try:
            mod = sm.tsa.statespace.SARIMAX( df['riders'],
                                            order=param,
                                            seasonal_order=param_seasonal,
                                            enforce_stationarity=False,
                                            enforce_invertibility=False
                                           )
            results = mod.fit()
            count += 1
            if count <= 5:
                print('ARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))
            
            if results.aic < select_candi:
                select_candi = results.aic
                param_candi = param
                param_seasonal_candi = param_seasonal
        except:
            continue
            
print(param_candi, param_seasonal_candi, select_candi)            

mod = sm.tsa.statespace.SARIMAX(
    df['riders'],
    order=(0, 1, 1),
    seasonal_order=(0, 1, 1, 12),
    enforce_stationarity=False,
    enforce_invertibility=False
)

results = mod.fit()

print( results.summary().tables[1] )

results.plot_diagnostics( figsize=(12, 10) )
plt.show()

pred = results.get_prediction(
    start=pd.to_datetime('1981-01-01'),
    dynamic=False
)

pred_ci = pred.conf_int()

# 관측 데이터  1973년 부터 끝까지 
ax = df['riders']['1973' : ].plot( label='observed', figsize=(12, 8) )

# 예측
pred.predicted_mean.plot(
    ax=ax,
    label='One-step ahead Forecast',
    alpha=.7
)

ax.fill_between(
    pred_ci.index,
    pred_ci['lower riders'],
    pred_ci['upper riders'],
    color='k',
    alpha=.2
)

ax.set_xlabel('Date')
ax.set_ylabel('Riders')
plt.legend(loc=5)

plt.show()

y_forecasted = pred.predicted_mean
y_truth = df['riders']['1981-01-01' : ]

# Compute the Mean Square Error
mse = (( y_forecasted - y_truth ) ** 2).mean()
print('The Mean Squared Error of our forecasts is {}'.format( round(mse, 2) ) ) 

pred_dynamic = results.get_prediction(
    start=pd.to_datetime('1981-01-01'),
    dynamic=True,
    full_results=True
)

pred_dynamic_ci = pred_dynamic.conf_int()

# 관측 데이터  1973년 부터 끝까지 
ax = df['riders']['1973' : ].plot( label='observed', figsize=(12, 8) )

# 예측
pred_dynamic.predicted_mean.plot(
    ax=ax,
    label='Dynamic Forcast'
)

ax.fill_between(
    pred_ci.index,
    pred_ci['lower riders'],
    pred_ci['upper riders'],
    color='k',
    alpha=.2
)

ax.fill_betweenx(
    ax.get_ylim(),
    pd.to_datetime('1981-01-01'),
    df.index[-1],
    alpha=.3
)

ax.set_xlabel('Date')
ax.set_ylabel('Riders')
plt.legend(loc=9)

plt.show()

# Extract the predicted and true values of our time series
y_forecasted = pred_dynamic.predicted_mean
y_truth = df['riders']['1981-01-01' : ]

# Compute the Mean Square Error
mse = (( y_forecasted - y_truth ) ** 2).mean()
print('The Mean Squared Error of our forecasts is {}'.format( round(mse, 2) ) ) 

# Get forecast some steps ahead in future
pred_uc = results.get_forecast( steps=20 )

# Get confidence intervals of forecasts
pred_ci = pred_uc.conf_int()

# 관측 데이터  1973년 부터 끝까지 
ax = df['riders']['1973' : ].plot( label='observed', figsize=(12, 8) )

# 예측
pred_uc.predicted_mean.plot(
    ax=ax,
    label='Forcast'
)

ax.fill_between(
    pred_ci.index,
    pred_ci['lower riders'],
    pred_ci['upper riders'],
    color='k',
    alpha=.2
)


ax.set_xlabel('Date')
ax.set_ylabel('Riders')
plt.legend(loc=9)

plt.show()

data = sm.datasets.co2.load_pandas()
y = data.data
y.head()

y = y['co2'].resample('MS').mean()
y = y.fillna(y.bfill())
y.head()

y.plot( figsize=(15, 6) )
plt.show()

# Define the p, d and q parameters to take any value between 0 and 2
p = d = q = range(0, 2)

# Generate all different combinations of p, q and d triplets
pdq = list( itertools.product(p, d, q))

# Generate all different combinations of seasonal p, q and d triplets
seasonal_pdq = [ (x[0], x[1], x[2], 12)  for x in pdq ]

print('Example of parameter combinations for Seasonal ARIMA ...')
print('SARIMAX:  {} x {}'.format(pdq[1], seasonal_pdq[1]) )
print('SARIMAX:  {} x {}'.format(pdq[1], seasonal_pdq[2]) )
print('SARIMAX:  {} x {}'.format(pdq[2], seasonal_pdq[3]) )
print('SARIMAX:  {} x {}'.format(pdq[2], seasonal_pdq[4]) )

warnings.filterwarnings("ignore")

select_candi = 10000000
param_candi = ( 0, 0, 0 )
param_seasonal_candi = ( 0, 0, 0)

count=0
end_count = len(pdq)

for param in pdq:   
    for param_seasonal in seasonal_pdq:
        try:
            mod = sm.tsa.statespace.SARIMAX(
                                           # 여기 !!! 
                                            y,
                                            order=param,
                                            seasonal_order=param_seasonal,
                                            enforce_stationarity=False,
                                            enforce_invertibility=False
                                           )
            results = mod.fit()
            count += 1
            if count <= 5:
                print('ARIMA{}x{}12 - AIC:{}'.format(param, param_seasonal, results.aic))
            
            if results.aic < select_candi:
                select_candi = results.aic
                param_candi = param
                param_seasonal_candi = param_seasonal
        except:
            continue
            
print(param_candi, param_seasonal_candi, select_candi)    

mod = sm.tsa.statespace.SARIMAX(
    y,
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 12),
    enforce_stationarity=False,
    enforce_invertibility=False
)

results = mod.fit()

print( results.summary().tables[1] )

results.plot_diagnostics( figsize=(12, 10) )
plt.show()

pred = results.get_prediction(
    start=pd.to_datetime('1998-01-01'),
    dynamic=False
)

pred_ci = pred.conf_int()

# 관측 데이터  1973년 부터 끝까지 
ax = y[ '1990' : ].plot( label='observed', figsize=(12, 8) )

# 예측
pred.predicted_mean.plot(
    ax=ax,
    label='One-step ahead Forecast',
    alpha=.7
)

ax.fill_between(
    pred_ci.index,
    pred_ci.iloc[ : , 0 ],
    pred_ci.iloc[ : , 1 ],
    color='k',
    alpha=.2
)

ax.set_xlabel('Date')
ax.set_ylabel('CO2 Levels')
plt.legend()

plt.show()

y_forecasted = pred.predicted_mean
y_truth = y['1998-01-01' : ]

# Compute the Mean Square Error
mse = (( y_forecasted - y_truth ) ** 2).mean()
print('The Mean Squared Error of our forecasts is {}'.format( round(mse, 2) ) ) 

pred_dynamic = results.get_prediction(
    start=pd.to_datetime('1998-01-01'),
    dynamic=True,
    full_results=True
)

pred_dynamic_ci = pred_dynamic.conf_int()

# 관측 데이터  1973년 부터 끝까지 
ax = y['1990' : ].plot( label='observed', figsize=(12, 8) )

# 예측
pred_dynamic.predicted_mean.plot(
    ax=ax,
    label='Dynamic Forcast'
)

ax.fill_between(
    pred_ci.index,
    pred_ci.iloc[ : , 0 ],
    pred_ci.iloc[ : , 1 ],
    color='k',
    alpha=.25
)

ax.fill_betweenx(
    ax.get_ylim(),
    pd.to_datetime('1998-01-01'),
    y.index[-1],
    alpha=.1,
    zorder=-1
)

ax.set_xlabel('Date')
ax.set_ylabel('CO2 Levels')
plt.legend()

plt.show()

# Extract the predicted and true values of our time series
y_forecasted = pred_dynamic.predicted_mean
y_truth = y['1998-01-01' : ]

# Compute the Mean Square Error
mse = (( y_forecasted - y_truth ) ** 2).mean()
print('The Mean Squared Error of our forecasts is {}'.format( round(mse, 2) ) ) 

# Get forecast some steps ahead in future
pred_uc = results.get_forecast( steps=250 )

# Get confidence intervals of forecasts
pred_ci = pred_uc.conf_int()

# 관측 데이터  1973년 부터 끝까지 
ax = y.plot( label='observed', figsize=(12, 8) )

# 예측
pred_uc.predicted_mean.plot(
    ax=ax,
    label='Forcast'
)

ax.fill_between(
    pred_ci.index,
    pred_ci.iloc[ : , 0 ],
    pred_ci.iloc[ : , 1 ],
    color='k',
    alpha=.25
)


ax.set_xlabel('Date')
ax.set_ylabel('CO2 Levels')
plt.legend()

plt.show()