For this week’s discussion, I picked the data set from Bank of America on yahoo finance: https://finance.yahoo.com/quote/BAC/history?p=BAC

#Reading our dataset

#Summary of our dataset

str(BOA)
## 'data.frame':    60 obs. of  7 variables:
##  $ Date     : chr  "2015-08-01" "2015-09-01" "2015-10-01" "2015-11-01" ...
##  $ Open     : num  17.9 15.9 15.5 16.9 17.5 ...
##  $ High     : num  18.1 16.5 17.4 18.1 17.9 ...
##  $ Low      : num  14.6 15.2 14.6 16.9 16.5 ...
##  $ Close    : num  16.4 15.6 16.8 17.4 16.8 ...
##  $ Adj.Close: num  15 14.2 15.4 16 15.4 ...
##  $ Volume   : num  2.33e+09 1.78e+09 1.85e+09 1.44e+09 1.77e+09 ...
head(BOA)
##         Date  Open  High   Low Close Adj.Close     Volume
## 1 2015-08-01 17.91 18.07 14.60 16.44  15.00615 2325559300
## 2 2015-09-01 15.95 16.48 15.25 15.58  14.22115 1777912500
## 3 2015-10-01 15.52 17.44 14.63 16.78  15.36581 1848594500
## 4 2015-11-01 16.90 18.09 16.87 17.43  15.96102 1439390500
## 5 2015-12-01 17.52 17.89 16.50 16.83  15.41159 1771496200
## 6 2016-01-01 16.45 16.59 12.94 14.14  12.98476 2648604300
tail(BOA)
##          Date  Open  High   Low Close Adj.Close     Volume
## 55 2020-02-01 33.00 35.45 27.70 28.50  28.12317 1072992600
## 56 2020-03-01 28.35 29.75 17.95 21.23  20.94929 2826499000
## 57 2020-04-01 19.93 25.32 19.51 24.05  23.88343 1636882500
## 58 2020-05-01 23.38 26.17 20.10 24.12  23.95295 1447546400
## 59 2020-06-01 24.28 29.01 23.02 23.75  23.58551 1801607100
## 60 2020-07-01 24.03 24.87 22.39 24.36  24.36000 1193407900

#Creating our time series

BOA.ts = ts(BOA$Adj.Close, frequency = 12, start = c(2015,8))
autoplot(BOA.ts) +xlab("Year") + ylab("Adjusting Closing Price in $")

#Splitting our time series into training and testing

train.ts = window(BOA.ts, end=c(2020,6))
test.ts = window(BOA.ts, start=c(2015,8))

#Building ARIMA Model and ETS models

M1.ARIMA = auto.arima(train.ts)
summary(M1.ARIMA)
## Series: train.ts 
## ARIMA(0,1,0) 
## 
## sigma^2 estimated as 4.329:  log likelihood=-124.79
## AIC=251.59   AICc=251.66   BIC=253.65
## 
## Training set error measures:
##                     ME    RMSE      MAE       MPE     MAPE      MASE
## Training set 0.1456674 2.06294 1.582365 0.3593268 7.018769 0.3124246
##                      ACF1
## Training set -0.009918847

#Building ETS model

M1.ETS = ets(train.ts, model = "ZZZ")
summary(M1.ETS)
## ETS(M,N,N) 
## 
## Call:
##  ets(y = train.ts, model = "ZZZ") 
## 
##   Smoothing parameters:
##     alpha = 0.9458 
## 
##   Initial states:
##     l = 14.848 
## 
##   sigma:  0.0909
## 
##      AIC     AICc      BIC 
## 326.8948 327.3312 333.1274 
## 
## Training set error measures:
##                     ME     RMSE      MAE      MPE     MAPE      MASE       ACF1
## Training set 0.1569204 2.065125 1.587724 0.397414 7.046455 0.3134828 0.03991386

#Forecasting ARIMA model

fcast.M1.ARIMA = forecast(M1.ARIMA,h=12)
plot(fcast.M1.ARIMA)
lines(test.ts, col="red")  
legend("topleft",lty=1,col=c("red","blue"),c("actual values","forecast"))

#Forecasting ETS Model

fcast.M1.ETS = forecast(M1.ETS,h=12)
plot(fcast.M1.ETS)
lines(test.ts, col="red")  
legend("topleft",lty=1,col=c("red","blue"),c("actual values","forecast"))

#Building GARCH model – For the first time

GBOA = garch(BOA.ts)
## 
##  ***** ESTIMATION WITH ANALYTICAL GRADIENT ***** 
## 
## 
##      I     INITIAL X(I)        D(I)
## 
##      1     3.486087e+01     1.000e+00
##      2     5.000000e-02     1.000e+00
##      3     5.000000e-02     1.000e+00
## 
##     IT   NF      F         RELDF    PRELDF    RELDX   STPPAR   D*STEP   NPRELDF
##      0    1  3.660e+02
##      1    2  2.146e+02  4.14e-01  5.48e+00  1.4e-02  2.0e+03  1.0e+00  5.50e+03
##      2    4  2.142e+02  1.79e-03  1.76e-03  7.2e-04  2.0e+00  5.0e-02  1.40e-01
##      3    5  2.136e+02  2.57e-03  2.23e-03  1.1e-03  2.2e+00  1.0e-01  2.56e-03
##      4    7  2.136e+02  6.60e-05  6.61e-05  5.9e-05  1.1e+01  5.3e-03  3.47e-04
##      5    9  2.136e+02  1.10e-04  1.10e-04  1.2e-04  2.1e+00  1.1e-02  2.66e-04
##      6   10  2.136e+02  1.31e-04  1.32e-04  2.4e-04  2.6e+00  2.1e-02  1.54e-04
##      7   12  2.136e+02  4.15e-06  4.15e-06  1.5e-05  1.1e+01  1.3e-03  2.15e-05
##      8   14  2.136e+02  6.94e-06  6.94e-06  3.1e-05  2.1e+00  2.6e-03  1.76e-05
##      9   16  2.136e+02  1.17e-06  1.17e-06  6.3e-06  1.8e+01  5.2e-04  1.07e-05
##     10   18  2.136e+02  2.13e-06  2.13e-06  1.3e-05  3.0e+00  1.0e-03  9.72e-06
##     11   20  2.136e+02  3.41e-06  3.41e-06  2.6e-05  2.0e+00  2.1e-03  7.64e-06
##     12   23  2.136e+02  5.70e-08  5.70e-08  5.4e-07  1.3e+02  4.2e-05  4.32e-06
##     13   25  2.136e+02  1.13e-07  1.13e-07  1.1e-06  1.8e+01  8.4e-05  4.47e-06
##     14   27  2.136e+02  2.23e-08  2.23e-08  2.2e-07  3.3e+02  1.7e-05  4.37e-06
##     15   29  2.136e+02  4.44e-08  4.44e-08  4.4e-07  4.2e+01  3.3e-05  4.36e-06
##     16   31  2.136e+02  8.81e-08  8.81e-08  8.8e-07  2.1e+01  6.7e-05  4.31e-06
##     17   34  2.136e+02  1.75e-09  1.75e-09  1.8e-08  4.0e+03  1.3e-06  4.23e-06
##     18   36  2.136e+02  3.50e-09  3.50e-09  3.5e-08  5.0e+02  2.7e-06  4.24e-06
##     19   38  2.136e+02  6.99e-09  6.99e-09  7.0e-08  2.5e+02  5.4e-06  4.24e-06
##     20   40  2.136e+02  1.40e-09  1.40e-09  1.4e-08  5.0e+03  1.1e-06  4.23e-06
##     21   42  2.136e+02  2.80e-10  2.80e-10  2.8e-09  2.5e+04  2.1e-07  4.23e-06
##     22   44  2.136e+02  5.59e-11  5.59e-11  5.6e-10  1.3e+05  4.3e-08  4.23e-06
##     23   46  2.136e+02  1.12e-11  1.12e-11  1.1e-10  6.3e+05  8.6e-09  4.23e-06
##     24   48  2.136e+02  2.24e-11  2.24e-11  2.2e-10  7.8e+04  1.7e-08  4.23e-06
##     25   50  2.136e+02  4.47e-12  4.47e-12  4.5e-11  1.6e+06  3.4e-09  4.23e-06
##     26   52  2.136e+02  8.94e-12  8.94e-12  9.0e-11  2.0e+05  6.8e-09  4.23e-06
##     27   54  2.136e+02  1.79e-11  1.79e-11  1.8e-10  9.8e+04  1.4e-08  4.23e-06
##     28   56  2.136e+02  3.58e-12  3.58e-12  3.6e-11  2.0e+06  2.7e-09  4.23e-06
##     29   58  2.136e+02  7.15e-13  7.15e-13  7.2e-12  9.8e+06  5.5e-10  4.23e-06
##     30   61  2.136e+02  1.46e-14  1.43e-14  1.4e-13  4.9e+08  1.1e-11  4.23e-06
##     31   63  2.136e+02  2.40e-15  2.86e-15  2.9e-14  2.4e+09  2.2e-12  4.23e-06
##     32   66  2.136e+02  2.28e-14  2.29e-14  2.3e-13  7.6e+07  1.8e-11  4.23e-06
##     33   69  2.136e+02  5.32e-16  4.58e-16  4.6e-15  1.5e+10  3.5e-13  4.23e-06
##     34   71  2.136e+02  9.32e-16  9.16e-16  9.2e-15  1.9e+09  7.0e-13  4.23e-06
##     35   72  2.136e+02 -4.68e+07  1.83e-15  1.8e-14  3.8e+09  1.4e-12  4.23e-06
## 
##  ***** FALSE CONVERGENCE *****
## 
##  FUNCTION     2.135736e+02   RELDX        1.840e-14
##  FUNC. EVALS      72         GRAD. EVALS      35
##  PRELDF       1.832e-15      NPRELDF      4.228e-06
## 
##      I      FINAL X(I)        D(I)          G(I)
## 
##      1    3.486220e+01     1.000e+00     1.401e-03
##      2    9.537147e-01     1.000e+00     1.127e-01
##      3    8.032921e-13     1.000e+00     2.551e-01
summary(GBOA)
## 
## Call:
## garch(x = BOA.ts)
## 
## Model:
## GARCH(1,1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## 0.7457 0.9471 0.9956 1.0618 1.2193 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)
## a0 3.486e+01   1.026e+03    0.034    0.973
## a1 9.537e-01   6.388e+00    0.149    0.881
## b1 8.033e-13   5.417e+00    0.000    1.000
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 1.8201, df = 2, p-value = 0.4025
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 0.18546, df = 1, p-value = 0.6667
plot(GBOA)

library(fGarch) #for GARCH model
## Loading required package: timeDate
## Loading required package: timeSeries
## 
## Attaching package: 'timeSeries'
## The following object is masked from 'package:zoo':
## 
##     time<-
## Loading required package: fBasics
library(rugarch)
## Loading required package: parallel
## 
## Attaching package: 'rugarch'
## The following object is masked from 'package:stats':
## 
##     sigma
library(quantmod)
## Loading required package: TTR
## 
## Attaching package: 'TTR'
## The following object is masked from 'package:fBasics':
## 
##     volatility
## Version 0.4-0 included new data defaults. See ?getSymbols.
library(rmgarch)
## 
## Attaching package: 'rmgarch'
## The following objects are masked from 'package:xts':
## 
##     first, last
## The following objects are masked from 'package:dplyr':
## 
##     first, last
M1.GARCH = garchFit(formula = ~garch(1,1), data = train.ts)
## 
## Series Initialization:
##  ARMA Model:                arma
##  Formula Mean:              ~ arma(0, 0)
##  GARCH Model:               garch
##  Formula Variance:          ~ garch(1, 1)
##  ARMA Order:                0 0
##  Max ARMA Order:            0
##  GARCH Order:               1 1
##  Max GARCH Order:           1
##  Maximum Order:             1
##  Conditional Dist:          norm
##  h.start:                   2
##  llh.start:                 1
##  Length of Series:          59
##  Recursion Init:            mci
##  Series Scale:              6.27541
## 
## Parameter Initialization:
##  Initial Parameters:          $params
##  Limits of Transformations:   $U, $V
##  Which Parameters are Fixed?  $includes
##  Parameter Matrix:
##                       U         V   params includes
##     mu     -37.02908711  37.02909 3.702909     TRUE
##     omega    0.00000100 100.00000 0.100000     TRUE
##     alpha1   0.00000001   1.00000 0.100000     TRUE
##     gamma1  -0.99999999   1.00000 0.100000    FALSE
##     beta1    0.00000001   1.00000 0.800000     TRUE
##     delta    0.00000000   2.00000 2.000000    FALSE
##     skew     0.10000000  10.00000 1.000000    FALSE
##     shape    1.00000000  10.00000 4.000000    FALSE
##  Index List of Parameters to be Optimized:
##     mu  omega alpha1  beta1 
##      1      2      3      5 
##  Persistence:                  0.9 
## 
## 
## --- START OF TRACE ---
## Selected Algorithm: nlminb 
## 
## R coded nlminb Solver: 
## 
##   0:     79.743710:  3.70291 0.100000 0.100000 0.800000
##   1:     77.303392:  3.81192 0.0800889 0.0947001 0.777029
##   2:     72.256050:  4.09168 0.0594709 0.121489 0.756411
##   3:     69.771178:  4.18718 0.0265619 0.130810 0.739854
##   4:     69.022880:  4.29806 0.0384650 0.153382 0.741550
##   5:     68.243546:  4.29522 0.0277395 0.154416 0.737282
##   6:     68.021533:  4.29738 0.0226070 0.156210 0.735014
##   7:     67.840062:  4.31641 0.0147524 0.167994 0.727283
##   8:     67.780420:  4.36378 0.0234580 0.173487 0.723136
##   9:     67.350435:  4.31263 0.0210322 0.180042 0.716275
##  10:     67.146102:  4.32431 0.0164987 0.192082 0.705725
##  11:     66.847463:  4.32239 0.0227089 0.204733 0.696350
##  12:     66.615504:  4.31998 0.0192022 0.216613 0.684821
##  13:     66.278567:  4.32056 0.0254987 0.241816 0.663091
##  14:     65.987920:  4.31958 0.0198269 0.261223 0.635921
##  15:     65.730118:  4.32874 0.0297725 0.272677 0.605741
##  16:     65.556880:  4.31724 0.0265267 0.305272 0.597714
##  17:     65.509834:  4.33676 0.0250803 0.306047 0.594648
##  18:     65.506776:  4.33826 0.0293260 0.308206 0.590523
##  19:     65.473951:  4.33668 0.0264988 0.308433 0.589205
##  20:     65.383790:  4.34157 0.0253648 0.319690 0.566664
##  21:     65.227224:  4.32624 0.0359681 0.337160 0.520657
##  22:     64.948348:  4.35401 0.0400496 0.424073 0.463584
##  23:     64.777615:  4.34120 0.0437692 0.490271 0.383110
##  24:     64.735483:  4.34122 0.0490963 0.536541 0.328511
##  25:     64.717027:  4.33182 0.0539522 0.571903 0.290392
##  26:     64.711803:  4.32491 0.0566985 0.585517 0.276344
##  27:     64.710991:  4.32295 0.0578448 0.588564 0.272757
##  28:     64.710917:  4.32290 0.0580753 0.588092 0.272847
##  29:     64.710911:  4.32303 0.0580771 0.587537 0.273197
##  30:     64.710911:  4.32306 0.0580653 0.587303 0.273317
##  31:     64.710910:  4.32304 0.0580635 0.587242 0.273314
##  32:     64.710910:  4.32303 0.0580659 0.587254 0.273288
##  33:     64.710910:  4.32302 0.0580668 0.587261 0.273281
## 
## Final Estimate of the Negative LLH:
##  LLH:  173.0726    norm LLH:  2.933434 
##         mu      omega     alpha1      beta1 
## 27.1287506  2.2867152  0.5872611  0.2732807 
## 
## R-optimhess Difference Approximated Hessian Matrix:
##                mu      omega      alpha1      beta1
## mu     -3.0948203 -0.5466344   0.8779323   2.159859
## omega  -0.5466344 -1.2524038  -1.8395149  -7.040958
## alpha1  0.8779323 -1.8395149 -50.0471588 -53.502921
## beta1   2.1598592 -7.0409581 -53.5029209 -94.621066
## attr(,"time")
## Time difference of 0.005975008 secs
## 
## --- END OF TRACE ---
## 
## 
## Time to Estimate Parameters:
##  Time difference of 0.04784489 secs
## Warning: Using formula(x) is deprecated when x is a character vector of length > 1.
##   Consider formula(paste(x, collapse = " ")) instead.
summary(M1.GARCH)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 1), data = train.ts) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x000000002d61f650>
##  [data = train.ts]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##       mu     omega    alpha1     beta1  
## 27.12875   2.28672   0.58726   0.27328  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu       27.1288      0.7210   37.627   <2e-16 ***
## omega     2.2867      1.7835    1.282   0.1998    
## alpha1    0.5873      0.3014    1.948   0.0514 .  
## beta1     0.2733      0.3044    0.898   0.3693    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -173.0726    normalized:  -2.933434 
## 
## Description:
##  Wed Jul 29 21:43:23 2020 by user: student 
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  4.481175  0.106396    
##  Shapiro-Wilk Test  R    W      0.924492  0.001298827 
##  Ljung-Box Test     R    Q(10)  87.95505  1.365574e-14
##  Ljung-Box Test     R    Q(15)  96.35712  6.37268e-14 
##  Ljung-Box Test     R    Q(20)  103.9441  2.463585e-13
##  Ljung-Box Test     R^2  Q(10)  6.05134   0.8109315   
##  Ljung-Box Test     R^2  Q(15)  23.84875  0.06770697  
##  Ljung-Box Test     R^2  Q(20)  26.08246  0.1631058   
##  LM Arch Test       R    TR^2   11.88778  0.4547341   
## 
## Information Criterion Statistics:
##      AIC      BIC      SIC     HQIC 
## 6.002461 6.143311 5.994023 6.057443

#Comparing the GARCH model’s results with the other models’ results

acc.ARIMA = accuracy(fcast.M1.ARIMA, test.ts)
acc.ETS = accuracy(fcast.M1.ETS, test.ts)

#Accuracy of ARIMA model

acc.ARIMA
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.1456674 2.062940 1.582365 0.3593268 7.018769 0.3124246
## Test set     0.7744870 0.774487 0.774487 3.1793389 3.179339 0.1529160
##                      ACF1
## Training set -0.009918847
## Test set               NA

#Accuracy of ETS model

acc.ETS
##                     ME      RMSE       MAE      MPE     MAPE      MASE
## Training set 0.1569204 2.0651246 1.5877243 0.397414 7.046455 0.3134828
## Test set     0.7551898 0.7551898 0.7551898 3.100122 3.100122 0.1491059
##                    ACF1
## Training set 0.03991386
## Test set             NA
#It appears that in terms of AIC and BIC the GARCH model does a way better job than the ETS or ARIMA model.
#Is it enough to say that the GARCH model is the best model ? Probably not because it would be nice to see how it does graphically and how it does in term of accuracy.
#With the current info, I would still rather use an ARIMA or ETS model over a GARCH model.

#Appendix - trying to forecast GARCH model

#M1.GARCH = garchFit(formula = ~garch(1,1),data=train.ts, cond.dist="QMLE")
#ug_spec = ugarchspec(mean.model = list(armaOrder = c(1,1)))
#if we want to specify a particular GARCH model, this is the code.
#ugfit = ugarchfit(spec = ug_spec, data = train.ts)
#ugfore = ugarchforecast(ugfit, n.ahead = 20)
#summary(ugfore)
#plot(ugfore)

`