在这里插入代码片# -*- coding: utf-8 -*-
"""
Created on Mon Oct 28 13:44:41 2019

@author: LIUJi
"""

import pandas as pd
import numpy as np
import pandas_datareader.data as web
import matplotlib.pyplot as plt
from arch import arch_model
import warnings
warnings.filterwarnings('ignore')
import statsmodels.api as sm
from matplotlib.pylab import rcParams
from statsmodels.tsa.arima_model import ARIMA

rcParams['figure.figsize'] = 16,5#定义图像大小
rcParams["font.size"]=20


Apple_stock_price=web.get_data_yahoo('AAPL')
closeprice = Apple_stock_price['Adj Close']
close_price = closeprice.reset_index(drop=True)

#Plot daily return 
returns = 100 * close_price.pct_change().dropna()
returns.plot()
plt.show()

#Compute the daily log return of the index
logreturn = np.log(close_price/close_price.shift(1))*100

#Check missing value for log return
missing_value = logreturn.isnull().sum()
print('Missing Valu3e Number =', missing_value)

#Fill missing value
logreturn_fill = logreturn.fillna(method = 'bfill')
log_return = logreturn_fill

#train_test_split
nlength=len(log_return)
ntrain=round(nlength*0.8)+1
ntest=nlength-ntrain
train_data=log_return[0:ntrain]
test_data=log_return[ntrain:]


#we use tra.diff()(differenced data), because this time series is unit root process.
fig,ax = plt.subplots(2,1,figsize=(20,10))
fig = sm.graphics.tsa.plot_acf(log_return, lags=50, ax=ax[0])
fig = sm.graphics.tsa.plot_pacf(log_return, lags=50, ax=ax[1])
plt.show()

#Fit SARIMA model
model_arima = ARIMA(train_data, order=(4,0,3))
fit = model_arima.fit()
res_df = fit.predict(start=0,end=1006)
plt.plot(res_df,'')
print(fit.summary())


#Fit GARCH model
model_garch = arch_model(train_data,mean = 'zero', vol='EGARCH', p=1, q=1, dist='t')
model_fit = model_garch.fit()
print(model_fit.summary())
model_fit.plot()


#Plot train data
plt.figure()
plt.plot(train_data,color="gray")
plt.title('Garch[1][1]')

#Get residenal and conditional_volatility
sigema = model_fit.conditional_volatility.values
sigema_df = pd.DataFrame(sigema)

#Caculate Mean, Lower, Upper
mean = res_df
lower_na = mean - 2*sigema_df
upper_na = mean + 2*sigema_df

lower = lower_na.iloc[:,:1]
upper = upper_na.iloc[:,:1]

plt.plot(mean, color = 'brown')
plt.plot(upper, color = 'green')
plt.plot(lower, color = 'red')

plt.show()


#Plot train data
plt.figure()
plt.plot(train_data[600:1007],color="gray")


#find variance
omega = model_fit.params[0]
alpha = model_fit.params[1]
beta = model_fit.params[2]

var=[]
log_array = np.array(log_return)
var.append(np.sqrt(omega + alpha*log_array[len(train_data)-1]**2 + beta*sigema[-1]**2))
for i in range(len(test_data)-1):
    temp = omega + alpha*log_array[len(train_data)-1]**2+ beta*var[i]**2
    var.append(np.sqrt(temp))

#Find upper and lower bound of var
var_df_upper = pd.DataFrame(var)
index_col = pd.DataFrame(np.arange(1007,1258))
var_df_upper['index_c'] = index_col
var_df_upper.set_index(['index_c'],inplace = True)
var_df_low = var_df_upper.apply(lambda x: x * -1)
var_df_low.columns = ['a']


#plot
x = np.linspace(1007,1257,251)
y_upper_na = np.array(var_df_upper.loc[x])
y_upper = y_upper_na.reshape(-1)
y_low_na = np.array(var_df_low.loc[x])
y_low = y_low_na.reshape(-1)
plt.plot(var_df_low,linewidth=3,color='blue',alpha=0.6)
plt.plot(var_df_upper,linewidth=3,color='blue',alpha=0.6)
plt.fill_between(x,y_upper,y_low,color='orange')

plt.title('Forecast of var')
plt.show()