Description

I used the Quantmod package to retrieve RCL’s daily prices from January 3, 2000 through March 23, 2020. The default source is Yahoo Finance and the data is structured as an xts object33 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. It also includes the initial impact from COVID-19.

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 build44 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%
• In 2020 the stock has dropped 74.8% YTD

Log Returns Distribution

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

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

Logged Returns
Autocorrelation Function
Partial Autocorrelation Function


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

Squared Logged Returns
Autocorrelation Function
Partial Autocorrelation Function


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.

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.750270

                      
Akaike       -4.750270
Bayes        -4.742561
Shibata      -4.750273
Hannan-Quinn -4.747570

I evaluated eight 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, and so I selected model six66 GARCH(1,1) with Skewed Student’s T Distribution. 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.02606268

sd(stdRCLret) 

[1] 1.016573

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, col="orange")

Standardized Residuals


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)

ACF of Standardized Residuals


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

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

ACF of Squared Standardized Residuals


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)

QQ Plot of Standardized Residuals


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)

Empirical Density of Standardized Residuals


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.000746    0.000253  2.944198 0.003238
omega   0.000003    0.000008  0.400919 0.688480
alpha1  0.051377    0.033446  1.536114 0.124510
beta1   0.947580    0.033250 28.498843 0.000000
skew    1.021471    0.017066 59.855060 0.000000
shape   4.420089    0.320901 13.774003 0.000000

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:77 However, the p-values for omega and alpha are >0.05. Since the coefficients are not significant they would likely be eliminated

\[R_t = 0.00074 + e_t\] \[e_t\sim N(0,\hat{\sigma}_t^2)\] \[\hat{\sigma}_t^2=0.000003+0.051377e^2_{t-1}+0.94758\hat{\sigma}^2_{t-1}\]

Model Predictions

#Unconditional Variance
garchuncvar <- uncvariance(garchfit6)
sqrt(garchuncvar) #long run SD 0.0546621

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

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

Estimated Volatilities


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

# Compare absolute returns versus predictions
plotvol <- plot(abs(RCLLogReturns), col = "steelblue1")
plotvol <- addSeries(garchvol, col = "black", on = 1, lwd=2)
plotvol

Estimated Volatility vs. Absolute Returns


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, 
                                  colorset="steelblue1",
                                  xlim=as.Date(c("2000-01-09", "2020-01-18")))

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

RCL annualized volatility, rolling 1 month


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).

Rolling window of 1 month (22 trading days)