In [1]:
import numpy as np
import matplotlib.pyplot as plt
N = 1000
SIGMA_EPS = 10
THETA = [0.9]
MU = 800
epsilon = np.random.normal(scale=SIGMA_EPS, size=(N,))
Y = np.ones((N,))*MU
for n in range(1,len(epsilon)):
Y[n] = MU + THETA[0]*epsilon[n-1] + epsilon[n]
plt.figure(figsize=(14,4))
plt.plot(Y)
plt.title('MA(1) model')
plt.grid()
plt.show()
In [10]:
from statsmodels.tsa.stattools import acf
# theoretical ACF
theoretical = np.zeros((100,))
theoretical[0] = 1
theoretical[1] = THETA[0]/(THETA[0]**2+1)
plt.plot(theoretical, 'o')
# from statsmodel
acf_sm, qstat, pval = acf(Y, nlags=100, qstat=True)
plt.plot(acf_sm, '-+')
plt.show()
In [11]:
from statsmodels.tsa.arima_model import ARIMA
model = ARIMA(Y, (0,0,1)).fit()
print (model.summary())
In [12]:
# Plot the residuals and test for their correlation
plt.figure(figsize=(12,4))
plt.plot(Y)
plt.plot(model.fittedvalues, '--+')
plt.plot(model.resid, '--')
plt.grid()
plt.legend(['Truth', 'Predicted', 'Residual'])
plt.show()
In [13]:
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 [14]:
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 [15]:
yhat, std_err, confint = model.forecast()
print ("Predicted value = {}, StdErr = {}, Confidence Interval = {}".format(yhat, std_err, confint))
In [ ]: