#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