1. Summary

2. Introduction

3. Explore the data

The dataset is downloaded from yahoo.com\(^{[6]}\). It contains 574 observations of weekly adjusted closing price for Apple Inc. from 2005 to 2015. They are the adjusted closing price of the last trading day of each week. Below is the plot of the data. The red lines show the mean of each plot.

Both the original price and log price show an increasing trend from 2005 to 2015. The returns appear to be stationary with a small positive mean 0.00546. As time increases, the price can eventually go to infinity with the cumulation of positive returns.
Since our goal is to analyze the volatility (variability of returns), we will focus on the demeaned returns in the following parts. Looking at the plot of demeaned returns, we can see high volatility at 2008, probably caused by the huge success of selling iPhone and financial crisis.

# 4. ARMA model: returns

Before we dig into the volatility, we apply a linear stationary model as a simple way to capture the dynamic of returns over time. First we plot the sample autocorrelation function of the demeaned returns.

The aurocorrelation reaches 1 at lag 0. For other positive lags, the ACF mainly fall within the two standard error limit dashed lines, suggesting the correlations are not significant at \(5\%\) level. This is similar to a white noise process.

Let \(Y_n\) denotes the demeaned return at time \(n\) and the ARMA(p,q) model is as follows:

\[ \begin{aligned} {\mathbb{E}}[Y_n] &= 0 \\ Y_n &= \sum_{i=1}^p\phi_i Y_{n-i} + \epsilon_n +\sum_{j=1}^q\psi_j\epsilon_{n-j} \\ {\epsilon_n} &\overset{iid}\sim N[0,\sigma^2] \end{aligned} \]

We construct the AIC table and select an ARMA(p,q) model with lowest AIC score.

MA0 MA1 MA2 MA3 MA4 MA5
AR0 -1838.60 -1836.61 -1834.62 -1833.87 -1832.01 -1830.78
AR1 -1836.61 -1834.61 -1832.62 -1833.20 -1831.22 -1829.93
AR2 -1834.62 -1832.93 -1837.59 -1831.17 -1829.20 -1846.74
AR3 -1833.99 -1833.04 -1839.92 -1830.34 -1827.84 -1845.39
AR4 -1832.13 -1831.04 -1829.06 -1836.00 -1832.26 -1836.34

The AIC table favors ARMA(0,0), since it is the simplest model with smallest AIC value. This suggests that the demeaned returns may behave like Gaussian white noise. We will check the normality of demeaned returns.

## 
##  Shapiro-Wilk normality test
## 
## data:  ret.de
## W = 0.97135, p-value = 3.796e-09

There are heavy tails in the Q-Q plot of demeaned returns. The p-value of shapiro-wilk test is very small, rejecting the null hypothesis that demeaned returns follow Gaussian distribution.

Now we do some diagnostics on residuals.

The residuals are basically not autocorrelated at \(5\%\) level. The loglikelihood of ARMA(0,0) is 921.3. The demeaned returns may behave like non-Gaussian white noise.

5. GARCH model: returns and volatility

Apart from returns in ARMA model, we also consider the volatility in GARCH model. We use GARCH model as a benchmark of POMP model. GARCH is a simpler and popular model of volatility clustering. Let \(Y_n\) denotes the demeaned return at time \(n\) and \(V_n\) denotes its conditional variance. The GARCH(p,q) model is as follows:

\[ \begin{aligned} Y_n &= \epsilon_n\sqrt{V_n} \\ V_n &= \alpha_0+\sum_{j=1}^p\alpha_j Y_{n-j}^2+\sum_{k=1}^q\beta_k V_{n-k} \\ \end{aligned} \] where \(\epsilon_n\) is a white noise process.

Aftering attempting several low-order GARCH models, we choose GARCH(1,1) based on coefficents significance and likelihood. Now we rewrite the model:

\[ \begin{aligned} Y_n &= \epsilon_n\sqrt{V_n} \\ V_n &= \alpha_0+\alpha_1 Y_{n-1}^2+\beta_1 V_{n-1} \\ \end{aligned} \]

A large \(Y_{n-1}^2\) or \(V_{n-1}\) results in a large value of \(V_n\). Thus a large \(Y_{n-1}\) tends to be followed by a large \(Y_n\) at the next time point. This corresponds to the volatility clustering\(^{[1]}\).

## 
## Call:
## garch(x = ret.de, order = c(1, 1), grad = "analytic", trace = FALSE)
## 
## Model:
## GARCH(1,1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.70251 -0.66627  0.05693  0.64990  3.31852 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)    
## a0 0.0001329   0.0000699    1.902   0.0572 .  
## a1 0.1167132   0.0223195    5.229  1.7e-07 ***
## b1 0.8325814   0.0457675   18.192  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 64.546, df = 2, p-value = 9.659e-15
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 0.023854, df = 1, p-value = 0.8773

The loglikelihood of GARCH(1,1) model is 938.189, which is larger than that of ARMA(0,0).

6. POMP model: returns, volatility and leverage

In previous section of GARCH model, we use a specific function to describe the volatility. GARCH model is useful for predicting volatility but its parameters are not explanatory. To study the correlation between return and volatility, we implement POMP model using stochastic equation to represent leverage.

We compare the maximum likelihood of two stochastic leverage models: fixed leverage and random-walk leverage. In the fixed leverage model, the leverge is constant over time. As an alternative, the latter model has time-varying leverage parameter following a random walk. Random walk is a simple idiosyncratic dynamic models with only one parameter of noise variance.

The random variation of leverage is idiosyncratic, i.e., independent of shocks to returns and of volatility. Ignoring random leverage could trigger both an excessive confidence in conclusions and biases. Random leverage could help to prevent poor hedging and control risk\(^{[7]}\).

First we build the POMP model folllowing the notations in Bretó (2014)\(^{[7]}\). Then we perform iterated filtering to estimate parameters maximizing log likelihood function. At each iteration, the particle filter is performed on a perturbed version of the model\(^{[9]}\). At last we use random starting values to approach MLE.

6.1 Fixed-leverage model

6.1.1 Build a model

State space variables: observable \(\{Y_n\}\), latent \(\{H_n\}\).
Parameters: \(\mu_h,\phi,\sigma_{\eta},\rho\).
\(\rho\): fixed leverage.
\(\{Y_n\}\): demeaned return on time \(n\), observable.
\(\{H_n\}\): log volatility, unobservable.
\(\{\epsilon_n\}\): return shock.

\[ \begin{align} Y_n &= \exp\{H_n/2\} \epsilon_n, \\ H_n &= \mu_h(1-\phi) + \phi H_{n-1} + \beta_{n-1}R_n\exp\{-H_{n-1}/2\} + \omega_n,\\ \end{align} \] where \[ \begin{align} \beta_n &= Y_n\sigma_\eta\sqrt{1-\phi^2} \\ \sigma_\omega &= \sigma_\eta\sqrt{1-\phi^2}\sqrt{1-\rho^2} \\ \epsilon_n &\overset{iid}\sim N[0,1] \\ \omega_n &\overset{iid}\sim N[0,\sigma_\omega^2] \\ \end{align} \]

6.1.2 Fit a model

The starting values of parameters and initial values of dynamic system combine repetitive experiments and the empirical results in Bretó (2014)\(^{[7]}\).

params_test_fix <- c(
     mu_h = -0.25,       
     phi = 0.98,     
     sigma_eta = 0.9,
    rho=-0.65,
      H_0=0
  )

Here is the summary of log likelihood using 5000 sequential Monte Carlo samples and 200 times of iterations. The values of log likelihood are smaller than those of ARMA and GARCH.

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  763.66  803.08  825.14  818.57  832.74  871.23

Below is the diagnostics plot of last iteration and the convergence diagnostics.

We can observe the latent variable volatility. After the initial decrease, it varies within a fixed range. However, the log likelihood does not converge.

The points are rather sparse in a neighbourhood of the likelihood surface. The results of likelihood inference are not satisfactoy in the fixed leverage model. This suggests that we should implement the random-walk leverage model with time-varying parameters.

6.2 Random-walk leverage model

6.2.1 Build a model

Simlar to 6.1.1, we build the random-walk leverage model. We use \(\{R_n\}\) instead of constant \(\rho\) to denote the leverage.
State space variables: observable \(\{Y_n\}\), latent \(\{G_n\},\{H_n\}\).
Parameters: \(\mu_h,\phi,\sigma_{\eta},\sigma_{\nu}\).
\(\{R_n\}\): leverage on time \(n\) as correlation between return on time \(n-1\) and the increase in the log volatility from time \(n-1\) to \(n\).It is a random walk on a transformed scale \([-1,1]\).
\(\{G_n\}\): usual Gaussian random walk leverage, unobservable.
\(\{Y_n\}\): demeaned return on time \(n\), observable.
\(\{H_n\}\): log volatility, unobservable.
\(\{\epsilon_n\}\): return shock.

\[ \begin{align} R_n &= \frac{\exp\{2G_n\} -1}{\exp\{2G_n\}+1} \\ Y_n &= \exp\{H_n/2\} \epsilon_n, \\ H_n &= \mu_h(1-\phi) + \phi H_{n-1} + \beta_{n-1}R_n\exp\{-H_{n-1}/2\} + \omega_n,\\ G_n &= G_{n-1}+\nu_n, \end{align} \] where \[ \begin{align} \beta_n &= Y_n\sigma_\eta\sqrt{1-\phi^2} \\ \sigma_\omega &= \sigma_\eta\sqrt{1-\phi^2}\sqrt{1-R_n^2} \\ \epsilon_n &\overset{iid}\sim N[0,1] \\ \nu_n &\overset{iid}\sim N[0,\sigma_{\nu}^2] \\ \omega_n &\overset{iid}\sim N[0,\sigma_\omega^2] \\ \end{align} \]

6.2.2 Fit a model

Here are the starting values of parameters and initial values of dynamic system.

params_test <- c(
     sigma_nu = 0.0086,  
     mu_h = -0.261,       
     phi = 0.98,     
     sigma_eta = 0.92,
     G_0 = 0,
     H_0=0
  )

The improvement lies in the values of log likelihood: they are larger than those in the fixed leverage model and close to the benchmark of GARCH. The time variation leverage outperforms fixed leverage.

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  890.00  894.29  902.14  902.31  910.03  914.61

As we can see, the likelihood function converges to a maximum value as iterations increase. The parameter \(\phi\) also converges to 1. The results in convergence also favors the random walk leverage model.

The maximum likelihood surface has denser points than the fixed leverage model. Smaller values of \(\sigma_{\nu}\) are more probable to maximize likelihood.

6.2.3 Likelihood Maximization using randomized starting values

Instead of the specific starting values in 6.2.2, we randomly select the starting values by uniform distribution within parameter vectors. If the conclusion of estimation are stable, this can support the finding of MLE via global search\(^{[2]}\). Here are the parameter values based on the parameter estimates and standard errors in Bretó (2014)\(^{[7]}\).

apple_box <- rbind(
sigma_nu=c(0.005,0.01),
 mu_h    =c(-0.4,0),
 phi = c(0.95,0.99),
 sigma_eta = c(0.8,1),
 G_0 = c(-1,1),
 H_0 = c(-0.5,0.5)
)

The best likelihood is slightly larger than that the best likelihood in 6.2.2, but the overall likelihood is smaller.

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  848.95  860.79  881.81  883.32  908.11  916.21

The likelihood converges faster in this case.

The points in the likelihood surface are more sparse than those of 6.2.2.

The diverse remote starting values can stably approach the maximum likelihood. This result has strengthened the MLE of random walk leverage model with specific starting values.

6.2.4 Profile likelihood

From the diagnostics plots given above, we can see that the parameter \(\phi\) can also converge to 1. This inspires us to investigate this parameter by constructing profile likelihood and rigorous confidence interval.

For each starting point, we find the maximum likelihood and its estimator. We fix \(\phi\) and let other parameters vary. We use the estimator of the largest likelihood value to construct the profile likelihood. Here is the expression of \(1-\alpha\) confidence interval of \(\phi\): \[ \{\phi:\max\{\ell^{\text{profile}}(\phi)\}-\ell^{\text{profile}}(\phi)< z_{\alpha/2} \} \]

##     lower     upper 
## 0.9500003 0.9603668

A \(95\%\) confidence interval of \(\phi\) is \([0.95,0.96]\). From the likelihood plot, there is no strong evidence showing that \(\phi\) is significant. Therefore, converging parameter does not necessarily contribute to converging likelihood.

7. Conclusion

8. Reference

[1] Ruey S. Tsay. (2002). Analysis of financial time series. Chicago, IL: John Wiley & Sons, Inc.
[2] http://ionides.github.io/531w16/notes15/notes15.html
[3] http://www.investopedia.com/terms/a/adjusted_closing_price.asp
[4] http://www.investopedia.com/terms/v/volatility.asp
[5] https://en.wikipedia.org/wiki/Apple_Inc.
[6] http://finance.yahoo.com/q/hp?s=AAPL&a=00&b=1&c=2005&d=11&e=31&f=2015&g=w
[7] Carles Bretó, 2014. On idiosyncratic stochasticity of financial leverage effects. Statistics and Probability Letters 91 (2014), 20–26.
[8] http://fic.wharton.upenn.edu/fic/papers/09/0905.pdf
[9] https://r-how.com/packages/pomp/mif