一、时间序列（2_2）AR、MA、ARMA简单编程（R语言）

library(TSA)
#eg.1
#ar（1）
set.seed(1)
series.all=arima.sim(n=48,list(ar=0.8))+100
future=window(series.all,start=41)
series=window(series.all,end=40)
model1=arima(series,order=c(1,0,0))
model1
result=plot(model1,n.ahead=8,ylab='Series&Forecasts',col=NULL,pch=19)
abline(h=coef(model1)[names(coef(model1))=='intercept'])
forecast=result$pred
cbind(future,forecast)
plot(model1,n.ahead=8,ylab='Seires,Forecasts, Actual&Limits',pch=19)
points(x=(41:48),y=future,pch=3)
abline(h=coef(model1)[names(coef(model1))=='intercept'])

#eg.2
#ma2

set.seed(2)
series.all=arima.sim(n=36,list(ma=c(-1,0.6)))+100
actual=window(series.all,start=33)
series=window(series.all,end=32)
model2=arima(series,order=c(0,0,2))
model2
result=plot(model2,n.ahead=4,ylab='Series&Forecasts',col=NULL,pch=19)
abline(h=coef(model2)[names(coef(model2))=='intercept'])
forecast=result$pred
cbind(actual,forecast)
plot(model2,n.ahead=4,ylab='Seires,Forecasts, Actual&Limits',type='o',pch=19)
points(x=(33:36),y=actual,pch=3)
abline(h=coef(model2)[names(coef(model2))=='intercept'])

例题三：analysis of U.S.GNP growth rate series
Notes:
Use 5% level in all tests.
The notation 𝜌i
is the lag } i autocorrelation coefficient

#### Problem 1:daily VIX index of the Chicago Board Options Exchange (CBOE)
#remove those missing values:
#unit root
#ACF and PACF of the growth rate

da=read.table("VIXCLS.txt",header=T)
head(da)
idx=c(1:nrow(da))[da$VALUE < 0.1]
VIX=da$VALUE[-idx]
lvix=log(VIX)

#a
par(mfcol=c(2,1))
tdx=c(1:length(VIX))/252+1990
#252 trading days a year
plot(tdx,VIX,xlab='year',ylab='VIX',type='l')
plot(tdx,lvix,xlab='year',ylab='ln(VIX)',type='l')

#b
require(fUnitRoots)
m1=ar(diff(lvix),method="mle")
m1$order
adfTest(lvix,lags=11,type="nc")
#no drift
#null hypothesis cannot be rejected.
#b1
library(tseries)
adf.test(lvix)
#c
t.test(diff(lvix))
#VIX should not have drift.Otherwise, there exists arbitrage opportunity.
#d
zt=diff(lvix)
par(mfcol=c(2,1))
acf(zt)
pacf(zt)
Box.test(zt,lag=10,type='Ljung')
#serial correlations in the zt series.

### Problem 2:zt of Problem 1.
#AR fit,forecast
#a
m2=arima(zt,order=c(11,0,0),include.mean=F)
m2
tsdiag(m2,gof=20)##对估计进行诊断,判断残差是否为白噪声 
m3=arima(lvix,order=c(11,1,0))
m3
p3=predict(m3,3)
p3
#c
up=p3$pred+1.96*p3$se
lb=p3$pred-1.96*p3$se
pintv=cbind(lb,up)
print(round(pintv,3))