Abstract

Periods of high or low volatility tend to persist. Stock market volatility is typically tied to fear. Fear of recession, fear that the markets are over-priced, fear of new government regulation, etc. These type of fears often result in a negative feedback loop. The stock market starts to drop, market experts predict further declines, people panic and sell, the losses result in more fear, people continue to sell, etc. The panic builds on itself until it is exhausted, often weeks or months later. It is for this reason that analyzing volatility is important.

Many models, like regression, assume constant variance. With assets returns this would not be optimal, as volatility varies during certain periods of time and depends on past variance. The GARCH process (generalized autoregressive conditional heteroskedasticity), being autoregressive, depends on past squared observations and past variances to model current variance.

In this report I will analyze several GARCH models, describe the selection process, and provide predictions for future observations.

Dataset Selection

Description

I used R’s Quantmod package to retrieve RCL’s daily prices from January 3, 2000 through November 29, 2019. The default source is Yahoo Finance and the data is structured as an xts object (extensible time series). I selected this date range because it includes two major recessions - the dotcom bubble, post 1999, and the “Great Recession” of 2008.

Motivation

I chose to analyze Royal Caribbean due to its extreme sensitivity to economic events. The cruise line industry is very capital intensive. The cost to build a ship can range from $250 million to over $1 billion and it can take several years to complete. For example, Royal Caribbean’s Oasis of the Seas cost $1.4 billion and took more than two years to build (the vessel was ordered in Feb 2006 and was completed in Oct 2009). Cruises are also a highly discretionary purchase. During recessions people often elect to skip vacations. This results in a steep decline in revenue, which is particularly troublesome for an industry with very high fixed costs. As a result the stock can be very volatile.

  • On 9/17/01 the stock dropped 50.62%. That was the first trading day since 9/11/01
  • In 2008 the stock dropped 70.4%
  • In 2009 the stock gained 101.1%

Data Exploration

library(PerformanceAnalytics)
library(quantmod)
library(rugarch)
library(astsa)
library(tseries)
library(knitr)
library(kableExtra)

#Retrieve RCL Data 
getSymbols("RCL",from = "2000-01-01",to = "2019-11-30")
[1] "RCL"
kable(as.data.frame(RCL[(nrow(RCL)-5):nrow(RCL),]))%>% kable_styling(full_width=F, position = "left")
RCL.Open RCL.High RCL.Low RCL.Close RCL.Volume RCL.Adjusted
2019-11-21 116.71 118.85 115.52 118.69 2252200 117.9644
2019-11-22 119.43 120.05 118.35 119.13 1642400 118.4017
2019-11-25 119.80 120.13 119.21 119.42 1375000 118.6899
2019-11-26 119.49 120.41 118.57 119.74 2379800 119.0080
2019-11-27 120.54 121.31 119.58 121.30 1218900 120.5585
2019-11-29 120.64 120.86 119.75 120.02 751600 119.2863 
# Plot adjusted closing price RCL
plot(Ad(RCL),col=4, main="RCL Adjusted Closing Price")

We can see from stock chart that the price is not stationary. There is an upward trend and the mean is not constant.

#Calculate and view RCL LogReturns
RCLLogReturns <- diff(log(Ad(RCL)))
RCLLogReturns <- RCLLogReturns[-1]
plot(RCLLogReturns, main="RCL Log Returns", col=6)

Instead we can model the daily returns. Since the daily change is typically small we can approximate with log returns.

\[\nabla(X_t)\approx r_t\]

Mean: 0.00026
Standard Deviation: 0.029

The mean is almost 0 and the daily standard deviation is almost 3%.

# Compute annualized standard deviation
sqrt(252) * sd(RCLLogReturns) #Annual SD 0.4660095
sqrt(252) * sd(RCLLogReturns["2001"]) #2001 annualized 0.7719925
sqrt(252) * sd(RCLLogReturns["2008"]) #2008 annualized 0.8202522
sqrt(252) * sd(RCLLogReturns["2013"]) #2008 annualized 0.2459386

However, the standard deviation differs over time. The annualized standard deviation is 46.6% (calculated as \(\sqrt{252}\) x daily \(\sigma\) , where \(252\) is the number of trading days in a year). The annualized standard deviation for 2001 was 77.2%, for 2008 it was 82%, and for 2013 it was 24.6%.

adf.test(RCLLogReturns) #ADF Test

    Augmented Dickey-Fuller Test

data:  RCLLogReturns
Dickey-Fuller = -16.036, Lag order = 17, p-value = 0.01
alternative hypothesis: stationary

Based on the Augmented Dickey-Fuller test, with a p-value <0.05, we can accept the alternative hypothesis that the log returns are stationary

Log Returns Distribution

#Chart RCL Returns distribution vs Normal Distribution
chart.Histogram(RCLLogReturns, methods = c("add.normal", "add.density"),
                colorset=c("gray","red","blue"),
                main="RCL Distribution in Red, Normal Distribution in Blue ")

The log returns are not normally distributed. They have a higher peak, fatter tails, and are negatively skewed. This suggests that using a skewed Student’s T distribution may be appropriate.

ACF and PACF review

acf2(RCLLogReturns) #ACF/PACF of logged returns

The ACF/PACF of the Log Returns do not show any clear patterns and would be considered white noise.

acf2(RCLLogReturns^2) #ACF/PACF of squared logged returns

The ACF/PACF of the squared log returns show many significant values indicating that there is a higher-order dependence structure. This indicates that a GARCH model would be appropriate.

Model Selection

Evaluation

#####################Constant mean, standard garch(1,1) model, sstd#############
garchspec6 <- ugarchspec(
  mean.model = list(armaOrder=c(0,0)),
  variance.model = list(model = "sGARCH",garchOrder=c(1,1)),
  distribution.model ="sstd")

garchfit6 <- ugarchfit(data=RCLLogReturns,
                       spec=garchspec6)

infocriteria(garchfit6) #AIC -4.763747
                      
Akaike       -4.763747
Bayes        -4.755939
Shibata      -4.763750
Hannan-Quinn -4.761011
Number ARMA order GARCH order Distribution AIC
1 (0,0) (1,1) Normal -4.611315
2 (0,0) (0,1) Normal -4.221188
3 (0,0) (1,0) Normal -4.370564
4 (0,0) (2,1) Normal -4.610882
5 (0,0) (2,2) Normal -4.614770
6 (0,0) (1,1) Skewed Student’s T -4.763747
7 (1,0) AR(1) (1,1) Skewed Student’s T -4.764838
8 (0,1) MA(1) (1,1) Skewed Student’s T -4.764867

I evaluated 8 models and used the AIC value to make my selection. Although the AR(1)+GARCH and MA(1)+GARCH had the lowest AIC values, their added complexity didn’t improve the model much from the more simple GARCH(1,1), and so I selected model #6. Increasing the GARCH order beyond (1,1) also did not improve the model much.

Residual Diagnostics

#Compute Standardized Residuals
stdRCLret <- residuals(garchfit6, standardize=TRUE)
mean(stdRCLret) 
[1] -0.02061841
sd(stdRCLret) 
[1] 1.018908

The mean of the standardized residuals is close to 0 and the standard deviation is close to 1. The standardized residuals are assumed to be i.i.d. N(0,1)

plot(stdRCLret, main="Standardized Residuals", col="orange")

There is one large outlier (when the stock dropped 50% in 2001). Other than that, the standardized residuals appear to have a constant mean and variance with no clear pattern.

#ACF of Standardized Residuals
plot(garchfit6, which=10)

The ACF of the standardized residuals would be considered a white noise pattern

#ACF of Squared Standardized Residuals
plot(garchfit6, which=11)

The ACF of the squared standardized residuals would be considered a white noise pattern. Although one is significant, the majority are not.

#QQ Plot standardized residuals
plot(garchfit6, which=9)

The QQ plot looks mostly normal, with the exception of a few outliers. The one large outlier dramatically alters the plot slope.

#Empirical Density of Standardized Residuals
plot(garchfit6, which=8)

The standardized residuals are not exactly normally distributed, but they are close (with the exception of a few outliers)

Model Coefficients

round(garchfit6@fit$matcoef,6) #Coefficents and P values
        Estimate  Std. Error    t value Pr(>|t|)
mu      0.000754    0.000280   2.695032 0.007038
omega   0.000003    0.000001   2.302800 0.021290
alpha1  0.042323    0.002821  15.002137 0.000000
beta1   0.955174    0.002085 458.074463 0.000000
skew    1.028420    0.019932  51.597626 0.000000
shape   4.589979    0.314693  14.585578 0.000000

Based on the p-values, which are all below 0.05, we can conclude that the model coefficients are all significant.

The GARCH(1,1) model with constant mean is represented by:

\[R_t =\mu +e_t\]

\[e_t\sim N(0,\sigma^2_t)\]

\[\sigma^2_t = \alpha_0 + \alpha_1e^2_{t-1}+\beta_1\sigma^2_{t-1}\] (Note: I have ignored the shape and skew for this analysis)

Based on the coefficients we arrive at the following model:

\[R_t = 0.000754 + e_t\] \[e_t\sim N(0,\hat{\sigma}_t^2)\] \[\hat{\sigma}_t^2=0.000003+0.042298e^2_{t-1}+0.955204\hat{\sigma}^2_{t-1}\] ##Model Predictions

The model predicts long-run standard deviation of 0.032(calculated as \(\sqrt{\frac{\alpha_0}{1-\alpha-\beta}}\), which is similar to the actual daily standard deviation of 0.029.

#Predicted Volatilities
garchvol <-sigma(garchfit6)
plot(garchvol, main="Estimated Volatilities",colorset=tim8equal)

This chart shows the predicted volatilities. We see the large spikes in 2001 and 2008 which correspond to the dramatic price declines.

# Compare absolute returns versus predictions
plotvol <- plot(abs(RCLLogReturns), col = "steelblue1",main="Estimated Volatility vs. Absolute Returns")
plotvol <- addSeries(garchvol, col = "black", on = 1, lwd=2)
plotvol

Here we can see the predicted volatilities relative to the absolute value of returns. The patterns match closely.

#Chart Rolling 1 month volatility versus predictions
plotvol<-chart.RollingPerformance(R = RCLLogReturns,
                                  width = 22,
                                  FUN = "sd.annualized",
                                  scale = 252, 
                                  main = "RCL annualized volatility, rolling 1 month",
                                  colorset="steelblue1",
                                  xlim=as.Date(c("2000-01-09", "2020-01-18")))

plotvol <- addSeries(sqrt(252)*garchvol, col = "red", on = 1, lwd=2)
plotvol

Rolling window of 1 month (22 trading days)

The model seems to fit the data very well, which can be visualized when the predictions are annualized. This chart compares the actual annualized volatility, calculated on a rolling 1 month basis (shown in Light Blue), compared to predicted annualized volatility (shown in Red).

Model Forecasting

For the next ten days the model predicts:

garchforecast6 <- ugarchforecast(fitORspec = garchfit6,n.ahead=1000)
cbind(head(fitted(garchforecast6),10),head(sigma(garchforecast6),10))
       2019-11-29 2019-11-29
T+1  0.0007544507 0.01534358
T+2  0.0007544507 0.01540885
T+3  0.0007544507 0.01547368
T+4  0.0007544507 0.01553808
T+5  0.0007544507 0.01560205
T+6  0.0007544507 0.01566560
T+7  0.0007544507 0.01572874
T+8  0.0007544507 0.01579147
T+9  0.0007544507 0.01585379
T+10 0.0007544507 0.01591572
# Plot annualized Forecasts, no simulation
fd <- rollapply(RCLLogReturns,width=22,FUN = "sd")
fd <-fd*sqrt(252)
fd=as.ts(fd)
ts.plot(fd,xlim=c(4000,6009),ylim=c(0,0.8), col=4, ylab="Annualized standard deviation",
        main="sigma forecast")
zs=sqrt(252)*sigma(garchforecast6)
af<-ts(zs, start=5009)
points(af,type="l",col=2)

This chart compares the actual annualized volatility (calculated on a rolling 1 month basis and shown in blue) compared to the annualized forecasted volatility (shown in red).

I have also used the model coefficients to forecast 4 iterations of 10 years of simulated returns, volatilities, and stock prices. (uses a random number generator to simulate the noise, \(\epsilon_t\))

#Simulate 10 years of returns
simgarchspec <- garchspec6
setfixed(simgarchspec) <-as.list(coef(garchfit6))
simgarch <- ugarchpath(spec=simgarchspec,m.sim=4,n.sim=10*252,rseed=12345)
simret<-fitted(simgarch)
plot.zoo(simret, main="10 years of Simulated Returns", col="purple")

plot.zoo(sigma(simgarch), main="10 years of Simulated Volatilities", col="springgreen4", lwd=1)

Future Price = \(P_{t+h}=P_te^{(r_{t+1}+r_{t+2}+...+r_{t+h})}\)

simprices<-exp(apply(simret,2,"cumsum"))
matplot(simprices,type="l",lwd=3, main="10 years of Simulated Prices")