HW2
Travis Compton
4/30/2021
Exercise 1: Seasonal ARIMA
Use the “Electricity.RData” examined in HW1 Exercise 2 (d)-(f). Fit the model using first order differenced log transformed series.
###(a) By visually checking, decide what SARIMA models seem appropriate, i.e., specify p, d, q and P,D, Q in SARIMA model, ARIMA(p, d, q) ∗ ARIMA(P,D, Q)𝑠. Choose the most appropriate two models and explain your answer.
We ended up choosing our two models because it seemed that there wasn’t much way to improve our ACF and Ljung box test.
log.electricity = log(electricity)
plot.ts(log.electricity)
diff.elec = diff(log.electricity, lag = 1)
plot.ts(diff.elec)
acf(diff.elec, lag = 60)
diff.elec2 = diff(diff.elec, lag = 12)
plot.ts(diff.elec2)
par(mfrow=c(1,2))
acf(diff.elec2, main = "acf");
pacf(diff.elec2, main = "pacf")
sarima(log.electricity, p=7, d=1, q=8, P=1, D=1, Q=1, S=12)## initial value -3.251628
## iter 2 value -3.461606
## iter 3 value -3.498319
## iter 4 value -3.553453
## iter 5 value -3.606488
## iter 6 value -3.629068
## iter 7 value -3.632268
## iter 8 value -3.633522
## iter 9 value -3.637607
## iter 10 value -3.638367
## iter 11 value -3.639661
## iter 12 value -3.643615
## iter 13 value -3.645674
## iter 14 value -3.646075
## iter 15 value -3.646822
## iter 16 value -3.647703
## iter 17 value -3.651195
## iter 18 value -3.653333
## iter 19 value -3.653786
## iter 20 value -3.654089
## iter 21 value -3.655301
## iter 22 value -3.655427
## iter 23 value -3.655538
## iter 24 value -3.655590
## iter 25 value -3.655676
## iter 26 value -3.655706
## iter 27 value -3.655721
## iter 28 value -3.655727
## iter 29 value -3.655734
## iter 30 value -3.655739
## iter 31 value -3.655744
## iter 32 value -3.655746
## iter 33 value -3.655749
## iter 34 value -3.655752
## iter 35 value -3.655754
## iter 36 value -3.655756
## iter 37 value -3.655756
## iter 38 value -3.655757
## iter 39 value -3.655758
## iter 40 value -3.655758
## iter 41 value -3.655758
## iter 42 value -3.655759
## iter 43 value -3.655759
## iter 44 value -3.655759
## iter 45 value -3.655759
## iter 46 value -3.655760
## iter 47 value -3.655760
## iter 48 value -3.655760
## iter 49 value -3.655761
## iter 50 value -3.655761
## iter 51 value -3.655761
## iter 52 value -3.655761
## iter 53 value -3.655761
## iter 54 value -3.655761
## iter 55 value -3.655761
## iter 55 value -3.655761
## iter 55 value -3.655761
## final value -3.655761
## converged
## initial value -3.647768
## iter 2 value -3.649519
## iter 3 value -3.651512
## iter 4 value -3.652541
## iter 5 value -3.652903
## iter 6 value -3.653371
## iter 7 value -3.653834
## iter 8 value -3.654343
## iter 9 value -3.655777
## iter 10 value -3.657573
## iter 11 value -3.659086
## iter 12 value -3.661817
## iter 13 value -3.662959
## iter 14 value -3.665449
## iter 15 value -3.668204
## iter 16 value -3.668882
## iter 17 value -3.669153
## iter 18 value -3.669790
## iter 19 value -3.670165
## iter 20 value -3.670421
## iter 21 value -3.670608
## iter 22 value -3.671094
## iter 23 value -3.671454
## iter 24 value -3.671881
## iter 25 value -3.671931
## iter 26 value -3.672165
## iter 27 value -3.673211
## iter 28 value -3.673493
## iter 29 value -3.673792
## iter 30 value -3.674121
## iter 31 value -3.674276
## iter 32 value -3.674489
## iter 33 value -3.674819
## iter 34 value -3.674848
## iter 35 value -3.675149
## iter 36 value -3.675336
## iter 37 value -3.675465
## iter 38 value -3.675503
## iter 39 value -3.675558
## iter 40 value -3.675598
## iter 41 value -3.675659
## iter 42 value -3.675682
## iter 43 value -3.675729
## iter 44 value -3.675759
## iter 45 value -3.675902
## iter 46 value -3.675948
## iter 47 value -3.676095
## iter 48 value -3.676214
## iter 49 value -3.676362
## iter 50 value -3.676487
## iter 51 value -3.676529
## iter 52 value -3.676593
## iter 53 value -3.676667
## iter 54 value -3.676740
## iter 55 value -3.676903
## iter 56 value -3.677061
## iter 57 value -3.677459
## iter 58 value -3.677541
## iter 59 value -3.677977
## iter 60 value -3.678613
## iter 61 value -3.679522
## iter 62 value -3.679675
## iter 63 value -3.680153
## iter 64 value -3.680524
## iter 65 value -3.680950
## iter 66 value -3.681104
## iter 67 value -3.681306
## iter 68 value -3.681362
## iter 69 value -3.681400
## iter 70 value -3.681423
## iter 71 value -3.681466
## iter 72 value -3.681481
## iter 73 value -3.681482
## iter 74 value -3.681484
## iter 75 value -3.681509
## iter 75 value -3.681461
## iter 76 value -3.681515
## iter 76 value -3.681498
## iter 76 value -3.681496
## final value -3.681515
## converged## Warning in sqrt(diag(fitit$var.coef)): NaNs produced
## Warning in sqrt(diag(fitit$var.coef)): NaNs produced
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), include.mean = !no.constant, transform.pars = trans, fixed = fixed,
## optim.control = list(trace = trc, REPORT = 1, reltol = tol))
##
## Coefficients:## Warning in sqrt(diag(x$var.coef)): NaNs produced## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ma1
## -0.1733 -0.0179 0.2833 -0.1141 -0.1448 -0.8956 -0.0733 -0.2369
## s.e. NaN NaN NaN NaN NaN NaN NaN NaN
## ma2 ma3 ma4 ma5 ma6 ma7 ma8 sar1
## -0.2927 -0.3647 0.1081 0.1608 0.9451 -0.3627 -0.2866 0.0773
## s.e. 0.0052 NaN 0.0082 NaN 0.0068 0.0064 NaN NaN
## sma1
## -0.991
## s.e. NaN
##
## sigma^2 estimated as 0.0005669: log likelihood = 866.55, aic = -1697.09
##
## $degrees_of_freedom
## [1] 366
##
## $ttable
## Estimate SE t.value p.value
## ar1 -0.1733 NaN NaN NaN
## ar2 -0.0179 NaN NaN NaN
## ar3 0.2833 NaN NaN NaN
## ar4 -0.1141 NaN NaN NaN
## ar5 -0.1448 NaN NaN NaN
## ar6 -0.8956 NaN NaN NaN
## ar7 -0.0733 NaN NaN NaN
## ma1 -0.2369 NaN NaN NaN
## ma2 -0.2927 0.0052 -55.9238 0
## ma3 -0.3647 NaN NaN NaN
## ma4 0.1081 0.0082 13.1359 0
## ma5 0.1608 NaN NaN NaN
## ma6 0.9451 0.0068 139.0317 0
## ma7 -0.3627 0.0064 -57.0873 0
## ma8 -0.2866 NaN NaN NaN
## sar1 0.0773 NaN NaN NaN
## sma1 -0.9910 NaN NaN NaN
##
## $AIC
## [1] -4.307342
##
## $AICc
## [1] -4.303211
##
## $BIC
## [1] -4.126975sarima(log.electricity, p=1, d=0, q=2, P=0, D=1, Q=1, S=12)## initial value -3.244613
## iter 2 value -3.311958
## iter 3 value -3.506558
## iter 4 value -3.512362
## iter 5 value -3.541675
## iter 6 value -3.553657
## iter 7 value -3.564956
## iter 8 value -3.570646
## iter 9 value -3.579411
## iter 10 value -3.589624
## iter 11 value -3.596638
## iter 12 value -3.603022
## iter 13 value -3.610762
## iter 14 value -3.614025
## iter 15 value -3.616584
## iter 16 value -3.616752
## iter 17 value -3.616912
## iter 18 value -3.617007
## iter 19 value -3.617022
## iter 20 value -3.617031
## iter 21 value -3.617031
## iter 22 value -3.617031
## iter 22 value -3.617031
## iter 22 value -3.617031
## final value -3.617031
## converged
## initial value -3.632154
## iter 2 value -3.634635
## iter 3 value -3.635718
## iter 4 value -3.636300
## iter 5 value -3.636440
## iter 6 value -3.636601
## iter 7 value -3.636727
## iter 8 value -3.636832
## iter 9 value -3.636925
## iter 10 value -3.636996
## iter 11 value -3.637030
## iter 12 value -3.637034
## iter 12 value -3.637034
## final value -3.637034
## converged
## $fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = constant, transform.pars = trans, fixed = fixed,
## optim.control = list(trace = trc, REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ma1 ma2 sma1 constant
## 0.9647 -0.3965 -0.2682 -0.8394 0.0021
## s.e. 0.0188 0.0545 0.0542 0.0307 0.0002
##
## sigma^2 estimated as 0.000668: log likelihood = 851.75, aic = -1691.5
##
## $degrees_of_freedom
## [1] 379
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.9647 0.0188 51.2963 0
## ma1 -0.3965 0.0545 -7.2749 0
## ma2 -0.2682 0.0542 -4.9467 0
## sma1 -0.8394 0.0307 -27.3519 0
## constant 0.0021 0.0002 10.6188 0
##
## $AIC
## [1] -4.282273
##
## $AICc
## [1] -4.281882
##
## $BIC
## [1] -4.222263(b)
Fit the two SARIMA models which you chose in part (a) using sarima function in R and interpret the result.
As mentioned in part (a) in either fit our lyung-box and acf plots for both of our fits look fairly poor as the p-values are very small and there are hardly any significant values on our ACF plot. The only improvement we saw was in our ljung box test as we increased our p and q values.
(c)
Check model diagnostics of two candidate models and choose the final one based on diagnostics plots and AIC/ BIC comparison
As it was the best model we were able to produce sarima(log.electricity, p=7, d=1, q=8, P=1, D=1, Q=1, S=12) is our final. Though it also had a higher AIC and BIC.
(d)
Make a prediction for next 4 values based on your final model. Make sure that you predict the values in original scale, not transformed values.
par(mfrow=c(1,1))
sarima.for(log.electricity, n.ahead=4, p=7, d=1, q=8, P=1, D=1, Q=1, S=12)
## $pred
## Jan Feb Mar Apr
## 2006 12.77689 12.65000 12.70064 12.61876
##
## $se
## Jan Feb Mar Apr
## 2006 0.02414625 0.02799875 0.02920825 0.03018942pred.elec = sarima.for(log.electricity, n.ahead=4, p=7, d=1, q=8, P=1, D=1, Q=1, S=12)
pred.elec## $pred
## Jan Feb Mar Apr
## 2006 12.77689 12.65000 12.70064 12.61876
##
## $se
## Jan Feb Mar Apr
## 2006 0.02414625 0.02799875 0.02920825 0.03018942exp(pred.elec$pred)## Jan Feb Mar Apr
## 2006 353941.9 311762.1 327956.2 302175.2Exercise 2: Exponential Smoothing Model
Use the data “gtemp” in astsa package analyzed in HW1 Exercise 4. Among 130 time series observations, we will use last 10 series (121-130) for the model validation and rest of them (1-120) for the model train. In other words, we aim to predict last 10 series based on fitted model using first 120 observations.
(a)
By visually checking first 120 time series plot, decide which ESM model seems the most appropriate. For example, answer should be like “Additive error + no trend + additive seasonality”. Explain your choice.
There does appear to be a trend so we can assert additive linear trend with little to no seasonality. So I would assume additive error + additive linear trend + additive seasonality.
data(gtemp)
train.gtemp = ts(gtemp[1:120])
test.gtemp = ts(gtemp[121:130])
plot.ts(train.gtemp)
(b)
Fit the model for the first 120 data which you chose in part (a). Use ets function with the specification of model input with three letters. Report autoplot() and interpret how your chosen model decompose the data.
My initial idea that seasonality existed within the time series plot was proved incorrect as as I began running different models, was given the error for Nonseasonal data so we could assume that we could have N for our seasonal classification. It appeared that additive error existed so we add that into our model the additive linear trend. Just to see we’ll also remove the linear trend to analyze the model. Lastly we add the Z's to our model to see if R will automatically select the best option. Looking at the summary of the model to automatically select the best option it chooses A,N,N.
gtemp.ets1 = ets(train.gtemp, model="ZZZ")
gtemp.ets2 = ets(train.gtemp, model="ANN")
gtemp.ets3 = ets(train.gtemp, model="AAN") summary(gtemp.ets1)## ETS(A,N,N)
##
## Call:
## ets(y = train.gtemp, model = "ZZZ")
##
## Smoothing parameters:
## alpha = 0.3775
##
## Initial states:
## l = -0.2644
##
## sigma: 0.0999
##
## AIC AICc BIC
## 25.59878 25.80568 33.96126
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.01414113 0.09905458 0.0810756 NaN Inf 0.9067666 0.1060173autoplot(gtemp.ets1)
(c)
Make a prediction for next 10 values based on your final model. Calculate the prediction error using test set. As we know the true Y from the test set, it is the mean of squared error, i.e., average of (Yt − 𝑌̂t)^2, for t=121,…,130.
When making a prediction for our final model we actually choose A,N,N and A,A,N to try and decide what is our best final model. Our test error for the A,N,N model produced a rate of 2.5% Our test error for the A,A,N model produced a rate of 1.9%
gtemp.ets2 = forecast(gtemp.ets2, h=10)
plot(gtemp.ets2)
mean((gtemp.ets2$mean - as.vector(test.gtemp))^2)## [1] 0.02519147gtemp.ets3 = forecast(gtemp.ets3, h=10)
plot(gtemp.ets3)
mean((gtemp.ets3$mean - as.vector(test.gtemp))^2)## [1] 0.01915633(d)
Now, you fit ARIMA for the first 120 data using the chosen model in Homework1 Exercise 4. Predict the next 10 values based on this ARIMA model. Just borrow your answer from HW 1 (i.e., no need to explain your model selection process), fit the model, and calculate the prediction error.
For the prediction error we return a value of 2.8%
case.fit = auto.arima(train.gtemp)
case.fit## Series: train.gtemp
## ARIMA(1,1,3) with drift
##
## Coefficients:
## ar1 ma1 ma2 ma3 drift
## -0.9200 0.4680 -0.6661 -0.3306 0.0053
## s.e. 0.0851 0.1157 0.0905 0.0928 0.0022
##
## sigma^2 estimated as 0.009165: log likelihood=112.44
## AIC=-212.88 AICc=-212.13 BIC=-196.21gtemp.arima = sarima(train.gtemp,1,1,3)## initial value -2.203453
## iter 2 value -2.321910
## iter 3 value -2.331302
## iter 4 value -2.333379
## iter 5 value -2.338872
## iter 6 value -2.342007
## iter 7 value -2.342347
## iter 8 value -2.342604
## iter 9 value -2.345418
## iter 10 value -2.349597
## iter 11 value -2.350915
## iter 12 value -2.353075
## iter 13 value -2.356011
## iter 14 value -2.357912
## iter 15 value -2.360197
## iter 16 value -2.360997
## iter 17 value -2.361160
## iter 18 value -2.361233
## iter 19 value -2.361250
## iter 20 value -2.361251
## iter 20 value -2.361251
## iter 20 value -2.361251
## final value -2.361251
## converged
## initial value -2.363785
## iter 2 value -2.363807
## iter 3 value -2.363812
## iter 4 value -2.363816
## iter 5 value -2.363823
## iter 6 value -2.363823
## iter 7 value -2.363825
## iter 8 value -2.363826
## iter 9 value -2.363827
## iter 10 value -2.363827
## iter 11 value -2.363827
## iter 11 value -2.363827
## iter 11 value -2.363827
## final value -2.363827
## converged
gtemp.pred = sarima.for(train.gtemp, n.ahead = 10, 1, 1, 3)
gtemp.pred## $pred
## Time Series:
## Start = 121
## End = 130
## Frequency = 1
## [1] 0.3006863 0.3384532 0.3561816 0.3500698 0.3658919 0.3615339 0.3757425
## [8] 0.3728690 0.3857118 0.3840950
##
## $se
## Time Series:
## Start = 121
## End = 130
## Frequency = 1
## [1] 0.09370193 0.10684668 0.11043014 0.11196768 0.11521256 0.11679882
## [7] 0.11977042 0.12139116 0.12413903 0.12578244mean((gtemp.pred$pred - as.vector(test.gtemp))^2)## [1] 0.02845975(e)
Compare prediction performance between ESM and ARIMA
ARIMA gives us an error of 2.8% ESM gives us an error of 1.9% at our lowest value.
Exercise 3: Exponential Smoothing Model
We repeat the same processes of Exercise 2 using the “Electricity.RData” analyzed in Exercise 1. Among 396 time series observations, we will use last 26 series (371-396) for the model validation and rest of them (1-370) for the model train. In other words, we aim to predict last 26 series based on fitted model using first 370 observations.
(a) By visually checking firt 370 time series plot, decide which ESM model seems the most appropriate. Explain your choice.
There does appear to be a trend so we can assert additive linear trend with apparent seasonality. So I would assume additive error + additive linear trend + additive seasonality.
train.elec = ts(electricity[1:370], freq=12)
test.elec = ts(electricity[371:396], freq=12)
plot.ts(train.elec)
(b)
Fit the model for the first 370 data which you chose in part (a). Use ets function with the specification of model input with three letters. Report autoplot() and interpret how your chosen model decompose the data.
Again, using Z’s to autoselect our model it chooses MAM as the best fit. To better compare different options we’ll use MMM and MAA against it. Looking at the plot MAM we see a clear seasonality trend. While level and slope are harder to interpret we can see a positive trend for level with a negative trend for slope.
elec.ets1 = ets(train.elec, model="ZZZ")
elec.ets2 = ets(train.elec, model="MAM")
elec.ets3 = ets(train.elec, model="MMM")
elec.ets4 = ets(train.elec, model="MAA")summary(elec.ets1)## ETS(M,Ad,M)
##
## Call:
## ets(y = train.elec, model = "ZZZ")
##
## Smoothing parameters:
## alpha = 0.6427
## beta = 5e-04
## gamma = 0.0067
## phi = 0.9786
##
## Initial states:
## l = 153445.5704
## b = 661.7965
## s = 1.0085 0.9265 0.9389 0.9852 1.1337 1.1436
## 1.0381 0.9614 0.9009 0.9693 0.935 1.0588
##
## sigma: 0.0276
##
## AIC AICc BIC
## 8678.796 8680.745 8749.239
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 624.7038 6691.059 4901.356 0.1873471 2.054026 0.5697099 0.07684258autoplot(elec.ets1)
(c)
Make a prediction for next 26 values based on your final model. Calculate the prediction error using validation set. As we know the true Y from the validation set, it is the mean of squared error, i.e., average of (Yt − 𝑌̂t)^2, for t=371,…,396.
Error rate for model 2: 1.42% Error rate for model 3: 1.58% Error rate for model 4: .90%
elec.pred2 = forecast(elec.ets2, h=26)
plot(elec.pred2)
mean((elec.pred2$mean - as.vector(test.elec))^2)/10^8## [1] 1.422288elec.pred3 = forecast(elec.ets3, h=26)
plot(elec.pred3)
mean((elec.pred3$mean - as.vector(test.elec))^2)/10^8## [1] 1.582313elec.pred4 = forecast(elec.ets4, h=26)
plot(elec.pred4)
mean((elec.pred4$mean - as.vector(test.elec))^2)/10^8## [1] 0.9036755(d)
Now, you fit SARIMA for the first 370 data using the chosen model in Exercise 1. Predict the next 26 values based on this SARIMA model. Just borrow your answer from Exercise 1, fit the model, and calculate the prediction error.
Again, to compare different models we’ll use the auto.arima function to pick the best combination and compare it against the one we chose in exercise 1. After computing the prediction error for the model we chose in exercise 1 we return a value of 1.17%. For the best chosen model we return a prediction error of .75%.
case.fit.elec = auto.arima(train.elec)
case.fit.elec## Series: train.elec
## ARIMA(1,0,2)(0,1,1)[12] with drift
##
## Coefficients:
## ar1 ma1 ma2 sma1 drift
## 0.9501 -0.4173 -0.2662 -0.6880 474.3421
## s.e. 0.0262 0.0597 0.0591 0.0352 61.6410
##
## sigma^2 estimated as 47261491: log likelihood=-3672.44
## AIC=7356.87 AICc=7357.11 BIC=7380.16elec.pred1 = sarima.for(train.elec, n.ahead = 26, p=7, d=1, q=8, P=1, D=1, Q=1, S=12)
elec.pred1## $pred
## Jan Feb Mar Apr May Jun Jul Aug
## 31
## 32 335021.8 290099.7 305061.8 292428.6 316725.2 341281.7 380116.4 385790.2
## 33 337143.3 294295.5 310294.6 296101.8 318793.1 345293.6 385876.2 389938.9
## Sep Oct Nov Dec
## 31 292052.0 318465.3
## 32 331747.6 308975.5 297001.7 324262.7
## 33 335035.2 314126.2 302347.9 327505.5
##
## $se
## Jan Feb Mar Apr May Jun Jul
## 31
## 32 7760.577 7979.583 8041.943 8159.018 8406.623 8552.441 8717.117
## 33 11021.795 11173.674 11296.690 11466.371 11674.998 11816.992 12027.884
## Aug Sep Oct Nov Dec
## 31 6579.846 7487.984
## 32 8977.935 9395.917 9628.958 10268.974 10694.055
## 33 12342.998 12685.762 12940.368 13673.042 14162.117elec.pred2 = sarima.for(train.elec, n.ahead = 26, p=1, d=0, q=2, P=0, D=1, Q=1, S=12)
elec.pred2## $pred
## Jan Feb Mar Apr May Jun Jul Aug
## 31
## 32 339757.0 299584.7 313093.9 295453.2 318173.1 343796.4 381626.5 384401.6
## 33 345904.7 305709.7 319197.2 301536.1 324236.5 349841.3 387653.8 390412.2
## Sep Oct Nov Dec
## 31 296090.8 323023.5
## 32 331395.8 312469.1 301741.5 329195.1
## 33 337390.5 318448.8 307706.8 335146.8
##
## $se
## Jan Feb Mar Apr May Jun Jul
## 31
## 32 7906.604 8058.420 8193.061 8312.740 8419.323 8514.398 8599.327
## 33 9666.075 9756.944 9838.258 9911.093 9976.391 10034.977 10087.574
## Aug Sep Oct Nov Dec
## 31 6826.466 7734.962
## 32 8675.286 8743.294 8804.239 9338.110 9564.413
## 33 10134.823 10177.291 10215.479 10637.927 10811.825mean((elec.pred1$pred - as.vector(test.elec))^2)/10^8## [1] 1.174907mean((elec.pred2$pred - as.vector(test.elec))^2)/10^8## [1] 0.7502739(e)
Compare prediction performance between ESM and ARIMA.
Best error rate for ESM: .90% Best error rate for Arima: .75%
Exercise 4: ARCH and GARCH Model
The “usdhkd.RData” contains daily USD/HKD (US dollar to Hong Kong dollar) exchange rate from January 1, 2005 to March 7, 2006
(a)
Generate the time series plot on exchange rate. Generate ACF and PACF from both raw data and squared data, respectively. Make comments based on your observation.
Looking at the time series plot it appears to be fairly stationary. Looking at the ACF we see a sharp spike at lag 1, and negative significant values at lag 2 and then tapers to zero. The PACF plot has a negative spike at lag 1 and 2 and tapers to zero after that. For the squared ACF plots we have the first two lags be significant a decay to zero and then a few smaller significant variables lag 7 and 8. For the PACF plot we have a spike at lag 1 with a the rest of the variables at zero with another small significant variable at lag 7, similar to what we see in the ACF plot.
plot.ts(usdhkd)
par(mfrow=c(1,2))
acf(usdhkd, main = "acf"); pacf(usdhkd, main = "pacf")
par(mfrow=c(1,2))
acf(usdhkd^2, main = "acf Squared"); pacf(usdhkd^2, main = "pacf Squared")
(b)
Fit total four mdoels – ARCH(1), ARCH(2), GARCH(1,1) and GARCH(2,2). Then choose the best one based on diagnostics test on residuals and AIC criteria. There is no one correct model and choose the one that you think the optimal
Most of the models looked very similar, but I ultimately ended up choosing GARCH(1,1) as the most optimal model. We come to this conclusion because it has the lowest AIC and BIC values.
ARCH(1). R appears to be correlated while R^2 appears to not be correlated. AIC: -4.52. BIC: -4.49
usd.arch1 = garchFit(~garch(1,0), data=usdhkd, trace=F)## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.summary(usd.arch1)##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(1, 0), data = usdhkd, trace = F)
##
## Mean and Variance Equation:
## data ~ garch(1, 0)
## <environment: 0x0000000025aa9fa0>
## [data = usdhkd]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1
## -0.00117994 0.00045482 0.49888267
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu -1.180e-03 1.053e-03 -1.120 0.262606
## omega 4.548e-04 4.296e-05 10.588 < 2e-16 ***
## alpha1 4.989e-01 1.336e-01 3.733 0.000189 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 977.8896 normalized: 2.268885
##
## Description:
## Sun May 02 15:45:37 2021 by user: FugBoi
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 4337.001 0
## Shapiro-Wilk Test R W 0.8287354 0
## Ljung-Box Test R Q(10) 20.34018 0.02619446
## Ljung-Box Test R Q(15) 26.32464 0.03474572
## Ljung-Box Test R Q(20) 30.61091 0.06054118
## Ljung-Box Test R^2 Q(10) 9.38141 0.4963288
## Ljung-Box Test R^2 Q(15) 10.35731 0.7966824
## Ljung-Box Test R^2 Q(20) 11.16525 0.9418239
## LM Arch Test R TR^2 9.46025 0.663197
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -4.523850 -4.495547 -4.523946 -4.512675ARCH(2). R appears to be correlated while R^2 appears to not be correlated. AIC: -4.60. BIC: -4.58
usd.arch2 = garchFit(~garch(2,0), data=usdhkd, trace=F)## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.summary(usd.arch2)##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(2, 0), data = usdhkd, trace = F)
##
## Mean and Variance Equation:
## data ~ garch(2, 0)
## <environment: 0x00000000230fceb8>
## [data = usdhkd]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 alpha2
## -0.00046702 0.00026598 0.54327969 0.48701087
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu -4.670e-04 8.699e-04 -0.537 0.591
## omega 2.660e-04 3.473e-05 7.658 1.89e-14 ***
## alpha1 5.433e-01 1.323e-01 4.108 3.99e-05 ***
## alpha2 4.870e-01 1.108e-01 4.395 1.11e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 997.0851 normalized: 2.313422
##
## Description:
## Sun May 02 15:45:37 2021 by user: FugBoi
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 1450.006 0
## Shapiro-Wilk Test R W 0.8744017 0
## Ljung-Box Test R Q(10) 18.6254 0.04528692
## Ljung-Box Test R Q(15) 26.25545 0.03542182
## Ljung-Box Test R Q(20) 28.93381 0.08907092
## Ljung-Box Test R^2 Q(10) 8.642008 0.5663814
## Ljung-Box Test R^2 Q(15) 11.90956 0.6858611
## Ljung-Box Test R^2 Q(20) 16.36083 0.6940018
## LM Arch Test R TR^2 9.24412 0.68195
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -4.608283 -4.570547 -4.608454 -4.593384GARCH(1,1). R appears to be correlated while R^2 appears to not be correlated. AIC: -4.61. BIC: -4.58
usd.garch1 = garchFit(~garch(1,1), data=usdhkd, trace=F)## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.summary(usd.garch1)##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(1, 1), data = usdhkd, trace = F)
##
## Mean and Variance Equation:
## data ~ garch(1, 1)
## <environment: 0x0000000029e7abd0>
## [data = usdhkd]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 beta1
## -1.6576e-04 1.5468e-05 2.4065e-01 8.0449e-01
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu -1.658e-04 8.902e-04 -0.186 0.8523
## omega 1.547e-05 1.904e-05 0.812 0.4167
## alpha1 2.406e-01 1.362e-01 1.767 0.0772 .
## beta1 8.045e-01 1.052e-01 7.648 2.04e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 999.3168 normalized: 2.3186
##
## Description:
## Sun May 02 15:45:37 2021 by user: FugBoi
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 2690.144 0
## Shapiro-Wilk Test R W 0.858535 0
## Ljung-Box Test R Q(10) 20.78275 0.02266047
## Ljung-Box Test R Q(15) 24.94071 0.05074532
## Ljung-Box Test R Q(20) 26.92712 0.1373269
## Ljung-Box Test R^2 Q(10) 3.768855 0.9571724
## Ljung-Box Test R^2 Q(15) 6.011632 0.9795517
## Ljung-Box Test R^2 Q(20) 8.392027 0.9889398
## LM Arch Test R TR^2 4.156909 0.9804356
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -4.618639 -4.580903 -4.618810 -4.603740GARCH(2,2). R appears to be correlated with the exception of two different test while R^2 appears to not be correlated. AIC: -4.61. BIC: -4.55
usd.garch2 = garchFit(~garch(2,2), data=usdhkd, trace=F)## Warning in sqrt(diag(fit$cvar)): NaNs produced## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
## Consider formula(paste(x, collapse = " ")) instead.summary(usd.garch2)##
## Title:
## GARCH Modelling
##
## Call:
## garchFit(formula = ~garch(2, 2), data = usdhkd, trace = F)
##
## Mean and Variance Equation:
## data ~ garch(2, 2)
## <environment: 0x000000002ad0d3e0>
## [data = usdhkd]
##
## Conditional Distribution:
## norm
##
## Coefficient(s):
## mu omega alpha1 alpha2 beta1 beta2
## -1.0068e-04 1.3247e-05 2.4447e-01 1.0000e-08 5.5623e-01 2.4114e-01
##
## Std. Errors:
## based on Hessian
##
## Error Analysis:
## Estimate Std. Error t value Pr(>|t|)
## mu -1.007e-04 8.768e-04 -0.115 0.90858
## omega 1.325e-05 9.320e-06 1.421 0.15519
## alpha1 2.445e-01 8.728e-02 2.801 0.00509 **
## alpha2 1.000e-08 NA NA NA
## beta1 5.562e-01 NA NA NA
## beta2 2.411e-01 NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log Likelihood:
## 1000.581 normalized: 2.321534
##
## Description:
## Sun May 02 15:45:37 2021 by user: FugBoi
##
##
## Standardised Residuals Tests:
## Statistic p-Value
## Jarque-Bera Test R Chi^2 2936.158 0
## Shapiro-Wilk Test R W 0.859941 0
## Ljung-Box Test R Q(10) 21.39029 0.01853056
## Ljung-Box Test R Q(15) 25.39566 0.04487562
## Ljung-Box Test R Q(20) 27.34351 0.1258733
## Ljung-Box Test R^2 Q(10) 3.013569 0.9811029
## Ljung-Box Test R^2 Q(15) 5.023137 0.9919287
## Ljung-Box Test R^2 Q(20) 7.216473 0.9959124
## LM Arch Test R TR^2 3.459938 0.9913337
##
## Information Criterion Statistics:
## AIC BIC SIC HQIC
## -4.615226 -4.558622 -4.615607 -4.592877(c)
Make a prediction for next 5 volatility based on your final model
(usd.pred = predict(usd.garch1, n.ahead=5))## meanForecast meanError standardDeviation
## 1 -0.0001657554 0.01384428 0.01384428
## 2 -0.0001657554 0.01468953 0.01468953
## 3 -0.0001657554 0.01552382 0.01552382
## 4 -0.0001657554 0.01635032 0.01635032
## 5 -0.0001657554 0.01717167 0.01717167predict(usd.garch1, n.ahead=5,plot=TRUE, conf=.9)
## meanForecast meanError standardDeviation lowerInterval upperInterval
## 1 -0.0001657554 0.01384428 0.01384428 -0.02293757 0.02260606
## 2 -0.0001657554 0.01468953 0.01468953 -0.02432788 0.02399637
## 3 -0.0001657554 0.01552382 0.01552382 -0.02570016 0.02536865
## 4 -0.0001657554 0.01635032 0.01635032 -0.02705965 0.02672813
## 5 -0.0001657554 0.01717167 0.01717167 -0.02841064 0.02807913