In [19]:
import numpy as np
import matplotlib.pyplot as plt
N = 1000
SIGMA_EPS = 10
THETA = [0.9, 0]
PHI = [0.7, 0]
MU = 10
epsilon = np.random.normal(scale=SIGMA_EPS, size=(N,))
Yp = np.ones((N,))*MU
for n in range(1,len(epsilon)):
Yp[n] = MU + PHI[0]*(Yp[n-1] - MU) + PHI[1]*(Yp[n-2] - MU) + THETA[0]*epsilon[n-1] + THETA[1]*epsilon[n-2] + epsilon[n]
plt.figure(figsize=(14,4))
plt.plot(Yp)
plt.title('ARMA(2,2) model')
plt.grid()
Y = Yp.cumsum()
plt.figure(figsize=(14,4))
plt.plot(Y)
plt.title('ARIMA(2,1,2) model')
plt.grid()
plt.show()
In [20]:
from statsmodels.tsa.stattools import acf
acf_sm, qstat, pval = acf(Yp, nlags=100, qstat=True)
plt.plot(acf_sm, '-+')
plt.show()
In [21]:
from statsmodels.tsa.arima_model import ARIMA
model = ARIMA(Yp, (1,0,1)).fit()
print (model.summary())
In [22]:
model = ARIMA(Y, (1,1,1)).fit()
print (model.summary())
In [24]:
# Plot the residuals and test for their correlation
plt.figure(figsize=(12,4))
plt.plot(Yp)
plt.plot(model.fittedvalues, '--+')
plt.plot(model.resid, '--')
plt.grid()
plt.legend(['Truth', 'Predicted', 'Residual'])
plt.show()
In [25]:
import pandas as pd
from statsmodels.api import qqplot
resid_df = pd.DataFrame(model.resid)
resid_df.plot(kind='kde')
qqplot(model.resid)
plt.show()
In [26]:
acf_res, qstat, pval = acf(model.resid, nlags=100, qstat=True)
plt.stem(acf_res,)
plt.hlines(0.05, 0,100, linestyle='dashed')
plt.hlines(-0.05, 0,100, linestyle='dashed')
plt.show()
In [27]:
yhat, std_err, confint = model.forecast()
print ("Predicted value = {}, StdErr = {}, Confidence Interval = {}".format(yhat, std_err, confint))
In [ ]: