#The data I used is monthly LEI collected by the US Conference Board. Leading Economic Index (LEI) contains ten components used to forecast the trend of US economic performance
#time series plot for the overall period
ts<-ts(LEI$LEI,start=c(1959,1),end=c(2021,2),frequency=12)
ts.test<-ts(LEI$LEI,start=c(2021,3),end=c(2022,2),frequency=12)
autoplot(ts)+
labs(title="US Leading Economic Index Time Series 1959-2021")

plot(decompose(ts,type=("additive")))

#the seasonal variation is constant over time
#checkresiduals for the training time series
checkresiduals(ts)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

#The residuals show an increasing trend over time
#All lags in ACF are statistical significant
#it doesn't show a normal distribution
#ETS(M,A,N)
mye<-ets(ts)
mye
## ETS(M,Ad,N)
##
## Call:
## ets(y = ts)
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 0.4741
## phi = 0.8438
##
## Initial states:
## l = 24.9394
## b = 0.1746
##
## sigma: 0.006
##
## AIC AICc BIC
## 3492.514 3492.627 3520.202
# AIC AICc BIC
#3492.514 3492.627 3520.202
plot(mye)

fc<-forecast(mye,h=12)
fc
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Mar 2021 111.5885 110.7349 112.4421 110.2831 112.8939
## Apr 2021 111.7475 110.2783 113.2167 109.5005 113.9945
## May 2021 111.8817 109.7926 113.9708 108.6867 115.0767
## Jun 2021 111.9949 109.2828 114.7071 107.8470 116.1428
## Jul 2021 112.0905 108.7580 115.4229 106.9940 117.1870
## Aug 2021 112.1711 108.2263 116.1159 106.1380 118.2041
## Sep 2021 112.2391 107.6934 116.7848 105.2871 119.1912
## Oct 2021 112.2965 107.1637 117.4293 104.4466 120.1464
## Nov 2021 112.3450 106.6403 118.0496 103.6205 121.0694
## Dec 2021 112.3858 106.1254 118.6463 102.8113 121.9604
## Jan 2022 112.4203 105.6203 119.2203 102.0206 122.8200
## Feb 2022 112.4494 105.1261 119.7728 101.2493 123.6495
ts.test%>%accuracy(fc$mean)
## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 85.99253 85.99276 85.99253 76.68778 76.68778 0.2154671 971.6032
# ME RMSE MAE MPE MAPE ACF1 Theil's U
#Test set 85.99253 85.99276 85.99253 76.68778 76.68778 0.2154671 971.6032
plot(fc)

#Differencing
adf.test(ts)
##
## Augmented Dickey-Fuller Test
##
## data: ts
## Dickey-Fuller = -3.7961, Lag order = 9, p-value = 0.01919
## alternative hypothesis: stationary
unitroot_ndiffs(ts)
## ndiffs
## 1
#needs 1 differencing
diff<-diff(ts,differences = 1)
autoplot(diff)+
labs(title="Differencing US LEI")

acf(diff)

pacf(diff)

autoplot(diff)

#The plot now is much stationary than before
#ARIMA (1,1,2)
#1 order of the autoregressive part
#1 degree of first differencing involved
#2 order of ,oving average part
arima<-auto.arima(ts)
arima
## Series: ts
## ARIMA(1,1,2) with drift
##
## Coefficients:
## ar1 ma1 ma2 drift
## 0.8884 -0.3655 -0.2351 0.1183
## s.e. 0.0319 0.0494 0.0414 0.0596
##
## sigma^2 estimated as 0.2109: log likelihood=-475.57
## AIC=961.15 AICc=961.23 BIC=984.22
#AIC=961.15 AICc=961.23 BIC=984.22
fc1<-forecast(arima,h=12)
fc1
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Mar 2021 111.5087 110.9202 112.0973 110.6087 112.4088
## Apr 2021 111.7663 110.6942 112.8385 110.1266 113.4061
## May 2021 112.0084 110.5208 113.4960 109.7333 114.2835
## Jun 2021 112.2366 110.3556 114.1176 109.3599 115.1134
## Jul 2021 112.4526 110.1898 114.7153 108.9920 115.9132
## Aug 2021 112.6577 110.0216 115.2937 108.6261 116.6892
## Sep 2021 112.8530 109.8510 115.8551 108.2618 117.4443
## Oct 2021 113.0398 109.6788 116.4008 107.8997 118.1800
## Nov 2021 113.2189 109.5062 116.9317 107.5407 118.8972
## Dec 2021 113.3913 109.3338 117.4487 107.1860 119.5966
## Jan 2022 113.5576 109.1627 117.9525 106.8362 120.2791
## Feb 2022 113.7186 108.9934 118.4437 106.4921 120.9451
ts.test%>%accuracy(fc1$mean)
## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 86.55913 86.56018 86.55913 76.80528 76.80528 0.8066062 425.3571
# ME RMSE MAE MPE MAPE ACF1 Theil's U
#Test set 86.55913 86.56018 86.55913 76.80528 76.80528 0.8066062 425.3571
plot(fc1)

#By comparing ETS with ARIMA, ARIMA tends to have a better prediction based on lower
#AICc, BIC, and the forecast plot is more similar to the original one, though it has a slightly higher RMSE
#GARCH
#armaOrder=1,1
mygarch1<-ugarchspec(variance.model = list(model="sGARCH",garchOrder=c(1,1)),
mean.model=list(armaOrder=c(1,1)),distribution.model = "std")
garch1<-ugarchfit(spec=mygarch1,data=ts)
garch1
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(1,0,1)
## Distribution : std
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 166.372624 8.147252 20.4207 0.000000
## ar1 0.998211 0.000277 3609.0785 0.000000
## ma1 0.308551 0.033902 9.1011 0.000000
## omega 0.025137 0.007629 3.2948 0.000985
## alpha1 0.610230 0.096532 6.3215 0.000000
## beta1 0.388770 0.051692 7.5208 0.000000
## shape 4.556097 0.592330 7.6918 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 166.372624 6.573531 25.3095 0.000000
## ar1 0.998211 0.000772 1293.4519 0.000000
## ma1 0.308551 0.038658 7.9815 0.000000
## omega 0.025137 0.023999 1.0474 0.294897
## alpha1 0.610230 0.116499 5.2380 0.000000
## beta1 0.388770 0.167350 2.3231 0.020174
## shape 4.556097 0.897144 5.0784 0.000000
##
## LogLikelihood : -340.0791
##
## Information Criteria
## ------------------------------------
##
## Akaike 0.93051
## Bayes 0.97381
## Shibata 0.93033
## Hannan-Quinn 0.94720
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 2.651 0.1035
## Lag[2*(p+q)+(p+q)-1][5] 25.169 0.0000
## Lag[4*(p+q)+(p+q)-1][9] 39.112 0.0000
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.003028 0.9561
## Lag[2*(p+q)+(p+q)-1][5] 0.008007 1.0000
## Lag[4*(p+q)+(p+q)-1][9] 0.013326 1.0000
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.002498 0.500 2.000 0.9601
## ARCH Lag[5] 0.006443 1.440 1.667 0.9997
## ARCH Lag[7] 0.009460 2.315 1.543 1.0000
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 8.261
## Individual Statistics:
## mu 2.1243
## ar1 1.3356
## ma1 0.1746
## omega 3.1353
## alpha1 1.2159
## beta1 4.2635
## shape 1.5255
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.69 1.9 2.35
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.5668 0.1176
## Negative Sign Bias 0.2834 0.7770
## Positive Sign Bias 0.1872 0.8516
## Joint Effect 2.6330 0.4517
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 15.93 0.6619
## 2 30 33.52 0.2574
## 3 40 43.81 0.2747
## 4 50 48.64 0.4877
##
##
## Elapsed time : 0.1747491
#Akaike 0.93051
#armaOrder=2,2 is the best with the lowest Akaike=0.498:
#look two months back to predict this monthly's LEI
mygarch2<-ugarchspec(variance.model = list(model="sGARCH",garchOrder=c(1,1)),
mean.model=list(armaOrder=c(2,2)),distribution.model = "std")
garch2<-ugarchfit(spec=mygarch2,data=ts)
garch2
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(2,0,2)
## Distribution : std
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 25.437785 0.209013 121.7041 0.000000
## ar1 1.906672 0.000413 4613.5319 0.000000
## ar2 -0.906362 0.000251 -3608.4060 0.000000
## ma1 -0.568977 0.038555 -14.7577 0.000000
## ma2 0.082556 0.038287 2.1563 0.031064
## omega 0.001692 0.000941 1.7987 0.072061
## alpha1 0.120970 0.025424 4.7581 0.000002
## beta1 0.878030 0.025422 34.5386 0.000000
## shape 4.899487 0.758010 6.4636 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 25.437785 0.042216 602.5628 0.000000
## ar1 1.906672 0.001750 1089.7175 0.000000
## ar2 -0.906362 0.000796 -1139.1427 0.000000
## ma1 -0.568977 0.044679 -12.7347 0.000000
## ma2 0.082556 0.043828 1.8836 0.059616
## omega 0.001692 0.000877 1.9296 0.053661
## alpha1 0.120970 0.026552 4.5560 0.000005
## beta1 0.878030 0.020951 41.9089 0.000000
## shape 4.899487 0.998372 4.9075 0.000001
##
## LogLikelihood : -176.7659
##
## Information Criteria
## ------------------------------------
##
## Akaike 0.49803
## Bayes 0.55371
## Shibata 0.49775
## Hannan-Quinn 0.51949
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.6381 0.424404
## Lag[2*(p+q)+(p+q)-1][11] 7.9652 0.001389
## Lag[4*(p+q)+(p+q)-1][19] 11.9922 0.193756
## d.o.f=4
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.07785 0.7802
## Lag[2*(p+q)+(p+q)-1][5] 0.32845 0.9807
## Lag[4*(p+q)+(p+q)-1][9] 0.38600 0.9993
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.001516 0.500 2.000 0.9689
## ARCH Lag[5] 0.023453 1.440 1.667 0.9983
## ARCH Lag[7] 0.035677 2.315 1.543 0.9999
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.4256
## Individual Statistics:
## mu 0.001165
## ar1 0.086176
## ar2 0.085220
## ma1 0.097237
## ma2 0.155731
## omega 0.639461
## alpha1 0.204677
## beta1 0.372448
## shape 0.119745
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 2.1 2.32 2.82
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.5094 0.6106
## Negative Sign Bias 0.3640 0.7159
## Positive Sign Bias 0.1337 0.8937
## Joint Effect 0.8641 0.8341
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 17.16 0.5788
## 2 30 22.42 0.8025
## 3 40 27.83 0.9087
## 4 50 40.46 0.8023
##
##
## Elapsed time : 0.3622899
plot(garch2,which='all')
##
## please wait...calculating quantiles...

fc3<-ugarchboot(garch2,n.ahead=10,method=c("Partial","Full")[1])
plot(fc3,which=2)

#it shows a gradually increasing trend forecast, though the actual values show a more steeper trend
#armaOrder=3,3
mygarch3<-ugarchspec(variance.model = list(model="sGARCH",garchOrder=c(1,1)),
mean.model=list(armaOrder=c(3,3)),distribution.model = "std")
garch3<-ugarchfit(spec=mygarch3,data=ts)
## Warning in arima(data, order = c(modelinc[2], 0, modelinc[3]), include.mean =
## modelinc[1], : possible convergence problem: optim gave code = 1
garch3
##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(3,0,3)
## Distribution : std
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 25.700962 0.172054 149.3775 0.000000
## ar1 0.981141 0.000349 2807.9132 0.000000
## ar2 0.851695 0.000330 2583.3617 0.000000
## ar3 -0.832214 0.000327 -2541.4235 0.000000
## ma1 0.357923 0.038521 9.2916 0.000000
## ma2 -0.443242 0.036672 -12.0867 0.000000
## ma3 0.096853 0.037743 2.5661 0.010284
## omega 0.001707 0.000944 1.8084 0.070550
## alpha1 0.123136 0.025868 4.7602 0.000002
## beta1 0.875864 0.025731 34.0396 0.000000
## shape 4.960388 0.769793 6.4438 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 25.700962 0.089515 287.1130 0.000000
## ar1 0.981141 0.000577 1701.4071 0.000000
## ar2 0.851695 0.000503 1692.1093 0.000000
## ar3 -0.832214 0.000493 -1686.3810 0.000000
## ma1 0.357923 0.045066 7.9423 0.000000
## ma2 -0.443242 0.043470 -10.1966 0.000000
## ma3 0.096853 0.043766 2.2130 0.026898
## omega 0.001707 0.000891 1.9170 0.055239
## alpha1 0.123136 0.027965 4.4033 0.000011
## beta1 0.875864 0.022770 38.4662 0.000000
## shape 4.960388 0.997116 4.9747 0.000001
##
## LogLikelihood : -177.4174
##
## Information Criteria
## ------------------------------------
##
## Akaike 0.50514
## Bayes 0.57319
## Shibata 0.50471
## Hannan-Quinn 0.53137
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.6719 0.412376
## Lag[2*(p+q)+(p+q)-1][17] 10.8619 0.001526
## Lag[4*(p+q)+(p+q)-1][29] 14.5756 0.533435
## d.o.f=6
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.08181 0.7749
## Lag[2*(p+q)+(p+q)-1][5] 0.33848 0.9795
## Lag[4*(p+q)+(p+q)-1][9] 0.39623 0.9993
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.001408 0.500 2.000 0.9701
## ARCH Lag[5] 0.023096 1.440 1.667 0.9984
## ARCH Lag[7] 0.034300 2.315 1.543 0.9999
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.5049
## Individual Statistics:
## mu 0.003772
## ar1 0.090120
## ar2 0.088834
## ar3 0.088200
## ma1 0.041535
## ma2 0.023847
## ma3 0.049724
## omega 0.633317
## alpha1 0.209410
## beta1 0.396772
## shape 0.093844
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 2.49 2.75 3.27
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.5285 0.5973
## Negative Sign Bias 0.3743 0.7083
## Positive Sign Bias 0.1227 0.9024
## Joint Effect 0.9098 0.8231
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 14.21 0.7710
## 2 30 13.81 0.9923
## 3 40 41.02 0.3819
## 4 50 54.00 0.2892
##
##
## Elapsed time : 0.5435951
#Akaike 0.50514