The purpose of this week’s discussion is to explore GARCH and fGARCH models using stock data. The company I’ve been analyzing off and on in these class discussions, YouGov PLC, should make a really interesting study in clustered volatility since they maintain a fairly high profile in election polling - which is defined in part by punctuated equilibrium.

The first step is to pull two years of share prices from Yahoo! Finance and create a ts() object.

library(quantmod)
library(forecast)
library(lubridate)
library(dplyr)
getSymbols('YOU.L', src = 'yahoo', from = Sys.Date() - years(2), to = Sys.Date())
[1] "YOU.L"
yougov <- ts(YOU.L, frequency = 365, start = as.numeric(Sys.Date() - years(2)))
yg.train <- window(yougov, end = 16743)
yg.test <- window(yougov, start = 16743)
head(yougov)
Time Series:
Start = c(16742, 1) 
End = c(16742, 6) 
Frequency = 365 
         YOU.L.Open YOU.L.High YOU.L.Low YOU.L.Close YOU.L.Volume YOU.L.Adjusted
16742.00    146.826    149.864   145.307      144.75        39025       142.9501
16742.00    147.838    151.889   147.838      148.00        21871       146.1596
16742.01    148.395    151.281   147.079      147.50        59111       145.6658
16742.01    150.370    150.370   146.826      145.00        11782       143.1969
16742.01    146.826    147.332   146.826      147.50       154268       145.6658
16742.01    147.311    150.370   146.826      148.50        64241       146.6534

We can get a general idea of what the error terms look like by plotting the additive decomposition of the data. There is a general upward trend in share prices and some seasonality

library(ggplot2)
library(yaztheme)
library(gridExtra)
grid.arrange(
  autoplot(yougov[,'YOU.L.Adjusted'])+
    theme_yaz() + 
    labs(title = 'Adjusted Closing Prices: YouGov PLC',
         y = 'Closing Price')+
    scale_x_continuous(breaks = c(16742,16743), labels = c('2 Yrs Ago',' 1 Yr Ago')),
  ggPacf(yougov[,'YOU.L.Adjusted'])+
    theme_yaz()+
    labs(title = 'PACF Plot'),
  nrow = 1
  )

Now we can fit a GARCH model to the adjusted closing price.The PACF plot indicates that a first-order autoregression component is appropriate. The prediction appears to be rougly an average of all points provided.

library(fGarch)
package <U+393C><U+3E31>fGarch<U+393C><U+3E32> was built under R version 3.4.2Loading required package: timeDate
package <U+393C><U+3E31>timeDate<U+393C><U+3E32> was built under R version 3.4.1Loading required package: timeSeries
package <U+393C><U+3E31>timeSeries<U+393C><U+3E32> was built under R version 3.4.2
Attaching package: <U+393C><U+3E31>timeSeries<U+393C><U+3E32>

The following object is masked from <U+393C><U+3E31>package:zoo<U+393C><U+3E32>:

    time<-

Loading required package: fBasics
package <U+393C><U+3E31>fBasics<U+393C><U+3E32> was built under R version 3.4.2

Rmetrics Package fBasics
Analysing Markets and calculating Basic Statistics
Copyright (C) 2005-2014 Rmetrics Association Zurich
Educational Software for Financial Engineering and Computational Science
Rmetrics is free software and comes with ABSOLUTELY NO WARRANTY.
https://www.rmetrics.org --- Mail to: info@rmetrics.org

Attaching package: <U+393C><U+3E31>fBasics<U+393C><U+3E32>

The following object is masked from <U+393C><U+3E31>package:TTR<U+393C><U+3E32>:

    volatility
garch.fit1 <- garchFit( ~ garch(1,0), data = yougov[,"YOU.L.Adjusted"])

Series Initialization:
 ARMA Model:                arma
 Formula Mean:              ~ arma(0, 0)
 GARCH Model:               garch
 Formula Variance:          ~ garch(1, 0)
 ARMA Order:                0 0
 Max ARMA Order:            0
 GARCH Order:               1 0
 Max GARCH Order:           1
 Maximum Order:             1
 Conditional Dist:          norm
 h.start:                   2
 llh.start:                 1
 Length of Series:          507
 Recursion Init:            mci
 Series Scale:              56.00741

Parameter Initialization:
 Initial Parameters:          $params
 Limits of Transformations:   $U, $V
 Which Parameters are Fixed?  $includes
 Parameter Matrix:
                      U         V   params includes
    mu     -38.54335183  38.54335 3.854335     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
    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 
     1      2      3 
 Persistence:                  0.1 


--- START OF TRACE ---
Selected Algorithm: nlminb 

R coded nlminb Solver: 

  0:     1118.3294:  3.85434 0.100000 0.100000
  1:     689.56978:  4.22767 0.717611 0.880496
  2:     681.35328:  4.17412 0.646898 0.904367
  3:     540.07806:  4.38805 0.0422732 0.918282
  4:     533.45342:  4.41607 0.0358957 0.918184
  5:     518.15965:  4.47211 0.0231452 0.918064
  6:     508.12203:  4.50405 0.0152774 0.918081
  7:     505.35835:  4.51317 0.0131053 0.918104
  8:     500.51250:  4.53216 0.00898666 0.918166
  9:     495.56437:  4.56464 0.00635638 0.918336
 10:     492.46822:  4.59265 0.0113416 0.918860
 11:     483.93064:  4.65381 0.00447871 0.922443
 12:     482.88980:  4.65609 0.00749640 0.922481
 13:     481.16584:  4.68787 0.00630607 0.925282
 14:     480.99507:  4.68798 0.00561837 0.925308
 15:     480.97691:  4.69008 0.00523483 0.925477
 16:     480.91012:  4.69429 0.00549967 0.926276
 17:     480.89001:  4.69805 0.00504348 0.927133
 18:     480.88702:  4.69903 0.00528592 0.927379
 19:     480.87870:  4.69756 0.00521663 0.927198
 20:     480.87270:  4.69684 0.00516878 0.927580
 21:     480.86853:  4.69919 0.00510935 0.928179
 22:     480.86630:  4.69896 0.00520273 0.928233
 23:     480.86358:  4.69872 0.00512805 0.928308
 24:     480.85861:  4.69854 0.00520986 0.928534
 25:     480.32966:  4.69308 0.00500850 0.982044
 26:     480.28422:  4.69514 0.00517604 0.984399
 27:     480.25935:  4.69699 0.00493222 0.986759
 28:     480.24246:  4.69855 0.00505710 0.989143
 29:     480.22835:  4.69961 0.00489173 0.991542
 30:     480.21244:  4.70177 0.00484211 0.998939
 31:     480.20931:  4.70057 0.00483926  1.00000
 32:     480.20918:  4.70082 0.00484080  1.00000
 33:     480.20918:  4.70083 0.00484056  1.00000

Final Estimate of the Negative LLH:
 LLH:  2521.13    norm LLH:  4.972642 
     mu   omega  alpha1 
263.281  15.184   1.000 

R-optimhess Difference Approximated Hessian Matrix:
               mu       omega       alpha1
mu     -1.1815006 -0.17918942    0.7537930
omega  -0.1791894 -0.09844452   -0.8388953
alpha1  0.7537930 -0.83889533 -212.3045407
attr(,"time")
Time difference of 0.01002598 secs

--- END OF TRACE ---


Time to Estimate Parameters:
 Time difference of 0.6368151 secs
garch_preds <- predict(garch.fit1, 60)%>%ts()
preds <- data.frame(yougov = yougov[,"YOU.L.Adjusted"],
                    garch = c(rep(NA,447),garch_preds[,'meanForecast']))%>%
  ts()
plot <- autoplot(preds)+
  labs(title = 'GARCH Predictions vs. Actual',
       y = 'Adjusted Closing Price')+
  theme_yaz()
plot

Comparing to other models from the class, my GARCH model performs poorly while other models are much more accurate.

ses.fit <- ses(yougov[,'YOU.L.Adjusted'])
ets.fit <- ets(yougov[,'YOU.L.Adjusted'])
I can't handle data with frequency greater than 24. Seasonality will be ignored. Try stlf() if you need seasonal forecasts.
arima.fit <- auto.arima(yougov[,'YOU.L.Adjusted'])
preds <- data.frame(yougov = yougov[,"YOU.L.Adjusted"],
                    garch = c(rep(NA,447),garch_preds[,'meanForecast']),
                    ses = c(rep(NA,447), ses.fit$fitted[448:507]),
                    ets = c(rep(NA,447), ets.fit$fitted[448:507]),
                    arima = c(rep(NA,447), arima.fit$fitted[448:507]))%>%
  ts()
autoplot(window(preds, 400, 507))+
  labs(y = 'Adjusted Closing Price',
       title = 'Models vs. Actual: YouGov Closing Prices')+
  theme_yaz()+
  scale_color_manual(values = yaz_cols)

LS0tDQp0aXRsZTogIlBSRURJQ1QgNDEzIC0gV2VlayA3IERpc2N1c3Npb24iDQphdXRob3I6ICdKb3NoIFlhem1hbicNCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KLS0tDQoNClRoZSBwdXJwb3NlIG9mIHRoaXMgd2VlaydzIGRpc2N1c3Npb24gaXMgdG8gZXhwbG9yZSBHQVJDSCBhbmQgZkdBUkNIIG1vZGVscyB1c2luZyBzdG9jayBkYXRhLiBUaGUgY29tcGFueSBJJ3ZlIGJlZW4gYW5hbHl6aW5nIG9mZiBhbmQgb24gaW4gdGhlc2UgY2xhc3MgZGlzY3Vzc2lvbnMsIFlvdUdvdiBQTEMsIHNob3VsZCBtYWtlIGEgcmVhbGx5IGludGVyZXN0aW5nIHN0dWR5IGluIGNsdXN0ZXJlZCB2b2xhdGlsaXR5IHNpbmNlIHRoZXkgbWFpbnRhaW4gYSBmYWlybHkgaGlnaCBwcm9maWxlIGluIGVsZWN0aW9uIHBvbGxpbmcgLSB3aGljaCBpcyBkZWZpbmVkIGluIHBhcnQgYnkgcHVuY3R1YXRlZCBlcXVpbGlicml1bS4gDQoNClRoZSBmaXJzdCBzdGVwIGlzIHRvIHB1bGwgdHdvIHllYXJzIG9mIHNoYXJlIHByaWNlcyBmcm9tIFlhaG9vISBGaW5hbmNlIGFuZCBjcmVhdGUgYSBgdHMoKWAgb2JqZWN0LiAgDQpgYGB7ciwgZWNobyA9IFRSVUUsIHdhcm5pbmc9RkFMU0UsIG1lc3NhZ2U9IEZBTFNFfQ0KbGlicmFyeShxdWFudG1vZCkNCmxpYnJhcnkoZm9yZWNhc3QpDQpsaWJyYXJ5KGx1YnJpZGF0ZSkNCmxpYnJhcnkoZHBseXIpDQpnZXRTeW1ib2xzKCdZT1UuTCcsIHNyYyA9ICd5YWhvbycsIGZyb20gPSBTeXMuRGF0ZSgpIC0geWVhcnMoMiksIHRvID0gU3lzLkRhdGUoKSkNCg0KeW91Z292IDwtIHRzKFlPVS5MLCBmcmVxdWVuY3kgPSAzNjUsIHN0YXJ0ID0gYXMubnVtZXJpYyhTeXMuRGF0ZSgpIC0geWVhcnMoMikpKQ0KeWcudHJhaW4gPC0gd2luZG93KHlvdWdvdiwgZW5kID0gMTY3NDMpDQp5Zy50ZXN0IDwtIHdpbmRvdyh5b3Vnb3YsIHN0YXJ0ID0gMTY3NDMpDQpoZWFkKHlvdWdvdikNCmBgYA0KDQpXZSBjYW4gZ2V0IGEgZ2VuZXJhbCBpZGVhIG9mIHdoYXQgdGhlIGVycm9yIHRlcm1zIGxvb2sgbGlrZSBieSBwbG90dGluZyB0aGUgYWRkaXRpdmUgZGVjb21wb3NpdGlvbiBvZiB0aGUgZGF0YS4gVGhlcmUgaXMgYSBnZW5lcmFsIHVwd2FyZCB0cmVuZCBpbiBzaGFyZSBwcmljZXMgYW5kIHNvbWUgc2Vhc29uYWxpdHkNCmBgYHtyLCBmaWcud2lkdGggPSA4LCBmaWcuaGVpZ2h0PTN9DQpsaWJyYXJ5KGdncGxvdDIpDQpsaWJyYXJ5KHlhenRoZW1lKQ0KbGlicmFyeShncmlkRXh0cmEpDQpncmlkLmFycmFuZ2UoDQogIGF1dG9wbG90KHlvdWdvdlssJ1lPVS5MLkFkanVzdGVkJ10pKw0KICAgIHRoZW1lX3lheigpICsgDQogICAgbGFicyh0aXRsZSA9ICdBZGp1c3RlZCBDbG9zaW5nIFByaWNlczogWW91R292IFBMQycsDQogICAgICAgICB5ID0gJ0Nsb3NpbmcgUHJpY2UnKSsNCiAgICBzY2FsZV94X2NvbnRpbnVvdXMoYnJlYWtzID0gYygxNjc0MiwxNjc0MyksIGxhYmVscyA9IGMoJzIgWXJzIEFnbycsJyAxIFlyIEFnbycpKSwNCiAgZ2dQYWNmKHlvdWdvdlssJ1lPVS5MLkFkanVzdGVkJ10pKw0KICAgIHRoZW1lX3lheigpKw0KICAgIGxhYnModGl0bGUgPSAnUEFDRiBQbG90JyksDQogIG5yb3cgPSAxDQogICkNCmBgYA0KDQpOb3cgd2UgY2FuIGZpdCBhIEdBUkNIIG1vZGVsIHRvIHRoZSBhZGp1c3RlZCBjbG9zaW5nIHByaWNlLlRoZSBQQUNGIHBsb3QgaW5kaWNhdGVzIHRoYXQgYSBmaXJzdC1vcmRlciBhdXRvcmVncmVzc2lvbiBjb21wb25lbnQgaXMgYXBwcm9wcmlhdGUuIFRoZSBwcmVkaWN0aW9uIGFwcGVhcnMgdG8gYmUgcm91Z2x5IGFuIGF2ZXJhZ2Ugb2YgYWxsIHBvaW50cyBwcm92aWRlZC4NCmBgYHtyfQ0KbGlicmFyeShmR2FyY2gpDQpnYXJjaC5maXQxIDwtIGdhcmNoRml0KCB+IGdhcmNoKDEsMCksIGRhdGEgPSB5b3Vnb3ZbLCJZT1UuTC5BZGp1c3RlZCJdKQ0KDQpnYXJjaF9wcmVkcyA8LSBwcmVkaWN0KGdhcmNoLmZpdDEsIDYwKSU+JXRzKCkNCnByZWRzIDwtIGRhdGEuZnJhbWUoeW91Z292ID0geW91Z292WywiWU9VLkwuQWRqdXN0ZWQiXSwNCiAgICAgICAgICAgICAgICAgICAgZ2FyY2ggPSBjKHJlcChOQSw0NDcpLGdhcmNoX3ByZWRzWywnbWVhbkZvcmVjYXN0J10pKSU+JQ0KICB0cygpDQpwbG90IDwtIGF1dG9wbG90KHByZWRzKSsNCiAgbGFicyh0aXRsZSA9ICdHQVJDSCBQcmVkaWN0aW9ucyB2cy4gQWN0dWFsJywNCiAgICAgICB5ID0gJ0FkanVzdGVkIENsb3NpbmcgUHJpY2UnKSsNCiAgdGhlbWVfeWF6KCkNCnBsb3QNCmBgYA0KDQpDb21wYXJpbmcgdG8gb3RoZXIgbW9kZWxzIGZyb20gdGhlIGNsYXNzLCBteSBHQVJDSCBtb2RlbCBwZXJmb3JtcyBwb29ybHkgd2hpbGUgb3RoZXIgbW9kZWxzIGFyZSBtdWNoIG1vcmUgYWNjdXJhdGUuIA0KYGBge3J9DQpzZXMuZml0IDwtIHNlcyh5b3Vnb3ZbLCdZT1UuTC5BZGp1c3RlZCddKQ0KZXRzLmZpdCA8LSBldHMoeW91Z292WywnWU9VLkwuQWRqdXN0ZWQnXSkNCmFyaW1hLmZpdCA8LSBhdXRvLmFyaW1hKHlvdWdvdlssJ1lPVS5MLkFkanVzdGVkJ10pDQoNCnByZWRzIDwtIGRhdGEuZnJhbWUoeW91Z292ID0geW91Z292WywiWU9VLkwuQWRqdXN0ZWQiXSwNCiAgICAgICAgICAgICAgICAgICAgZ2FyY2ggPSBjKHJlcChOQSw0NDcpLGdhcmNoX3ByZWRzWywnbWVhbkZvcmVjYXN0J10pLA0KICAgICAgICAgICAgICAgICAgICBzZXMgPSBjKHJlcChOQSw0NDcpLCBzZXMuZml0JGZpdHRlZFs0NDg6NTA3XSksDQogICAgICAgICAgICAgICAgICAgIGV0cyA9IGMocmVwKE5BLDQ0NyksIGV0cy5maXQkZml0dGVkWzQ0ODo1MDddKSwNCiAgICAgICAgICAgICAgICAgICAgYXJpbWEgPSBjKHJlcChOQSw0NDcpLCBhcmltYS5maXQkZml0dGVkWzQ0ODo1MDddKSklPiUNCiAgdHMoKQ0KYXV0b3Bsb3Qod2luZG93KHByZWRzLCA0MDAsIDUwNykpKw0KICBsYWJzKHkgPSAnQWRqdXN0ZWQgQ2xvc2luZyBQcmljZScsDQogICAgICAgdGl0bGUgPSAnTW9kZWxzIHZzLiBBY3R1YWw6IFlvdUdvdiBDbG9zaW5nIFByaWNlcycpKw0KICB0aGVtZV95YXooKSsNCiAgc2NhbGVfY29sb3JfbWFudWFsKHZhbHVlcyA9IHlhel9jb2xzKQ0KYGBgDQoNCg==