library(xts)
library(zoo)
library(tseries)
library(timeSeries)
library(forecast)
library(quantmod)
library(fGarch)
library(ggplot2)
setwd("C:/Users/User/Google 雲端硬碟/政治大學ECO/時間序列/學號對應資料序列")
hw<-read.csv("106258017.csv")
hwx<-as.ts(hw$x)
hwy<-as.ts(hw$y)
x值
(a).
下圖為時間序列圖，呈現上升的趨勢。
簡單統基分析中，Mean為18.769，Q1為14.493，Q2為18.245，Q3為23.483，最小值為8.493，最大值為28.894。偏態係數為0.08029649，呈現右偏態。峰態係數為-1.194696，呈現低闊峰。
plot(hwx)
summary(hw$x)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8.493 14.439 18.245 18.769 23.482 28.894
skewness(hw$x)
## [1] 0.08029649
## attr(,"method")
## [1] "moment"
kurtosis(hw$x)
## [1] -1.194696
## attr(,"method")
## [1] "excess"
(b).
下圖為ACF和PACF圖，從PACF圖中，t期與t-1期有很大的cov，因此初步估計為ARMA(1,0)模型。
(c).
從ACF和PACF中並沒有看出明顯的季節性因素，從時間序列圖中也沒有明顯的結構性改變。
(d).
在ARMA(1,0)中，係數為統計顯著，然而在Ljung-Box test，我們發現殘差仍然有序列相關。
由於PACF中，落後第三期有顯著相關，因此，我再估了ARMA(3,0)，係數仍然為統計顯著，殘差不具序列相關，同時AIC和BIC值都比ARMA(1,0)較優。因此選擇ARMA(3,0)。
model1<-arima(x = hwx,order = c(1,0,0))
summary(model1)
##
## Call:
## arima(x = hwx, order = c(1, 0, 0))
##
## Coefficients:
## ar1 intercept
## 0.9890 18.6217
## s.e. 0.0081 3.4139
##
## sigma^2 estimated as 0.6805: log likelihood = -369.85, aic = 745.7
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.04429824 0.8249181 0.6254039 0.03927822 3.598412 0.9996485
## ACF1
## Training set 0.1168566
r1<-model1$residuals
Box.test(x = r1,lag = 4,type = "Ljung-Box")
##
## Box-Ljung test
##
## data: r1
## X-squared = 16.342, df = 4, p-value = 0.002593
Box.test(x = r1,lag = 8,type = "Ljung-Box")
##
## Box-Ljung test
##
## data: r1
## X-squared = 25.287, df = 8, p-value = 0.00139
Box.test(x = r1,lag = 12,type = "Ljung-Box")
##
## Box-Ljung test
##
## data: r1
## X-squared = 29, df = 12, p-value = 0.00394
AIC(model1)
## [1] 745.6989
BIC(model1)
## [1] 756.8103
model2<-arima(x = hwx,order = c(3,0,0))
summary(model2)
##
## Call:
## arima(x = hwx, order = c(3, 0, 0))
##
## Coefficients:
## ar1 ar2 ar3 intercept
## 1.1367 -0.3392 0.1926 18.6790
## s.e. 0.0567 0.0847 0.0568 3.5848
##
## sigma^2 estimated as 0.6449: log likelihood = -361.86, aic = 733.72
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.04560523 0.803038 0.6072341 0.05802163 3.478929 0.9706058
## ACF1
## Training set 0.0175752
r2<-model2$residuals
Box.test(x = r2,lag = 4,type = "Ljung-Box")
##
## Box-Ljung test
##
## data: r2
## X-squared = 3.0705, df = 4, p-value = 0.5461
Box.test(x = r2,lag = 8,type = "Ljung-Box")
##
## Box-Ljung test
##
## data: r2
## X-squared = 9.7518, df = 8, p-value = 0.2829
Box.test(x = r2,lag = 12,type = "Ljung-Box")
##
## Box-Ljung test
##
## data: r2
## X-squared = 12, df = 12, p-value = 0.4456
AIC(model2)
## [1] 733.7241
BIC(model2)
## [1] 752.243
(e).
透過LM test我們發現殘差平方和其落後項具有統計上的顯著相關，具有ARCH effect，並且也從ACF和PACF中發現此現象。
r2_sq<-r2^2
lm.test<-lm(r2_sq[5:300]~r2_sq[4:299]+r2_sq[3:298]
+r2_sq[2:297]+r2_sq[1:296])
summary(lm.test)
##
## Call:
## lm(formula = r2_sq[5:300] ~ r2_sq[4:299] + r2_sq[3:298] + r2_sq[2:297] +
## r2_sq[1:296])
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.2442 -0.4918 -0.2761 0.1496 9.2827
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.42414 0.08804 4.818 2.34e-06 ***
## r2_sq[4:299] 0.28287 0.05867 4.821 2.30e-06 ***
## r2_sq[3:298] 0.16999 0.06054 2.808 0.00533 **
## r2_sq[2:297] -0.12832 0.06053 -2.120 0.03485 *
## r2_sq[1:296] 0.02339 0.05866 0.399 0.69034
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.18 on 291 degrees of freedom
## Multiple R-squared: 0.1256, Adjusted R-squared: 0.1136
## F-statistic: 10.45 on 4 and 291 DF, p-value: 6.389e-08
acf(r2_sq)
pacf(r2_sq)
(f).
透過殘差平方的ACF和PACF圖，我們對殘差估計GARCH(0,2)，其殘差不具序列關。
藉由此模型，將mean equation的殘差標準化，發現其平方項與自身的落後項不具統計顯著，因此暫不考慮Leverage effects。
modelgarch<-garchFit(data = r2,formula = ~garch(2,0))
summary(modelgarch)
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(2, 0), data = r2)
##
## Mean and Variance Equation:
## data ~ garch(2, 0)
## <environment: 0x0000000024bbfb90>
## [data = r2]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 alpha2
## 0.022636 0.434556 0.219405 0.088911
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.02264 0.04326 0.523 0.6008
## omega 0.43456 0.05513 7.883 3.11e-15 ***
## alpha1 0.21941 0.09300 2.359 0.0183 *
## alpha2 0.08891 0.06771 1.313 0.1892
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## -346.8443 normalized: -1.156148
##
## Description:
## Sat May 19 00:15:38 2018 by user: User
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 8.115506 0.01728782
## Shapiro-Wilk Test R W 0.9922924 0.1221714
## Ljung-Box Test R Q(10) 10.73062 0.3788864
## Ljung-Box Test R Q(15) 23.38672 0.07626709
## Ljung-Box Test R Q(20) 30.45523 0.06280556
## Ljung-Box Test R^2 Q(10) 11.01206 0.3565792
## Ljung-Box Test R^2 Q(15) 12.79923 0.6178019
## Ljung-Box Test R^2 Q(20) 41.90307 0.002847144
## LM Arch Test R TR^2 8.947378 0.7074184
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## 2.338962 2.388346 2.338613 2.358726
ht<-modelgarch@h.t
s2<-r2/sqrt(ht)
le.test<-lm(I(s2^2)[5:300]~s2[4:299]+s2[3:298]+
s2[2:297]+s2[1:296])
summary(le.test)
##
## Call:
## lm(formula = I(s2^2)[5:300] ~ s2[4:299] + s2[3:298] + s2[2:297] +
## s2[1:296])
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5369 -0.8794 -0.5625 0.1813 13.1426
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.004880 0.095490 10.523 <2e-16 ***
## s2[4:299] -0.004408 0.095731 -0.046 0.963
## s2[3:298] 0.135534 0.095313 1.422 0.156
## s2[2:297] -0.047520 0.095420 -0.498 0.619
## s2[1:296] -0.046356 0.095488 -0.485 0.628
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.635 on 291 degrees of freedom
## Multiple R-squared: 0.008453, Adjusted R-squared: -0.005176
## F-statistic: 0.6202 on 4 and 291 DF, p-value: 0.6484
(g).
最終我們建立ARMA(3,0)+GARCH(0,2)模型。
我們發現mean equation係數中，比較原本ARMA和後來加上GARCH的模型：
ARMA(3,0)
ar1 ar2 ar3
1.1367 -0.3392 0.1926
ARMA(3,0)+GARCH(0,2)
ar1 ar2 ar3
1.000000 -0.130922 0.105190
我們可以很明顯的發現雖然正負號沒變，但數值差異很大，同時在Ljung-Box test，我們發現殘差具有序列相關。
model3<-garchFit(data = hwx,formula = ~arma(3,0)+garch(2,0))
summary(model3)
##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~arma(3, 0) + garch(2, 0), data = hwx)
##
## Mean and Variance Equation:
## data ~ arma(3, 0) + garch(2, 0)
## <environment: 0x00000000223e1b98>
## [data = hwx]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu ar1 ar2 ar3 omega alpha1
## 0.508645 1.000000 -0.130922 0.105190 0.377737 0.344461
## alpha2
## 0.082171
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.50864 0.16695 3.047 0.00231 **
## ar1 1.00000 0.06320 15.822 < 2e-16 ***
## ar2 -0.13092 0.08458 -1.548 0.12166
## ar3 0.10519 0.05714 1.841 0.06564 .
## omega 0.37774 0.05491 6.879 6.04e-12 ***
## alpha1 0.34446 0.11652 2.956 0.00311 **
## alpha2 0.08217 0.06496 1.265 0.20588
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## -344.2129 normalized: -1.147376
##
## Description:
## Sat May 19 00:15:38 2018 by user: User
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 3.569202 0.167864
## Shapiro-Wilk Test R W 0.9949972 0.4405192
## Ljung-Box Test R Q(10) 12.74653 0.238191
## Ljung-Box Test R Q(15) 26.3499 0.03450183
## Ljung-Box Test R Q(20) 34.03145 0.02591273
## Ljung-Box Test R^2 Q(10) 12.79341 0.2354531
## Ljung-Box Test R^2 Q(15) 16.07316 0.377206
## Ljung-Box Test R^2 Q(20) 37.80216 0.009361472
## LM Arch Test R TR^2 11.58273 0.4797454
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## 2.341419 2.427841 2.340363 2.376005
(i).
下圖為最終估計的 estimated conditional variance 序列。
plot(as.ts(model3@h.t))
(j).
下圖為實際值、one-step-ahead forecast 預測值與預測誤差的比較。
f.v<-data.frame(time=201:300)
f.v$predict<-NA
i<-1
for(i in 1:100){
model.forecast<-garchFit(
data = hwx[1:(199+i)],
formula = ~arma(2,1)+garch(2,0))
pre.v<-predict(object = model.forecast,n.ahead=1)
f.v[i,2]<-pre.v$meanForecast
}
f.v$real<-hwx[201:300]
f.v$error<-f.v$real-f.v$predict
ggplot(data = f.v)+
geom_line(aes(x=f.v$time,y=f.v$predict),color="red")+
geom_line(aes(x=f.v$time,y=f.v$real))+
geom_line(aes(x=f.v$time,y=f.v$error),color="blue")+
labs(title="實際值(黑線)、預測值(紅線)、預測誤差(藍線",x="期數",y="數值")
(k).
下圖為實際值、no-change forecast 預測值與預測誤差的比較。
f.v2<-data.frame(time=201:300)
f.v2$predict<-hwx[200:299]
f.v2$real<-hwx[201:300]
f.v2$error<-f.v2$real-f.v2$predict
ggplot(data = f.v2)+
geom_line(aes(x=f.v2$time,y=f.v2$predict),color="red")+
geom_line(aes(x=f.v2$time,y=f.v2$real))+
geom_line(aes(x=f.v2$time,y=f.v2$error),color="blue")+
labs(title="實際值(黑線)、預測值(紅線)、預測誤差(藍線",x="期數",y="數值")
(l).
根據Diebold-Mariano Test，縱然看起來透過ARMA(3,0)+GARCH(0,2)的MSPE比較小，但統計上我們仍然無法拒絕“H0: 兩者預測一樣好”，前者預測方式在統計上並無顯著優於後者。
dm.test(f.v2$error,f.v$error,alternative = "less")
##
## Diebold-Mariano Test
##
## data: f.v2$errorf.v$error
## DM = -0.42885, Forecast horizon = 1, Loss function power = 2,
## p-value = 0.3345
## alternative hypothesis: less