\documentclass[notitlepage]{article} \usepackage{../../Math545} \usepackage{listings} \usepackage{numprint} \usepackage{enumerate} \usepackage{bbm} \usepackage{bm} \usepackage{verbatim} \usepackage{amsfonts,amsmath,amssymb} \newtheorem{alg}{Algorithm} \newtheorem{ex}{Example} \newtheorem{note}{Note} \newenvironment{proof}{\par\noindent{\normalfont{\bf{Proof }}}} {\hfill \small$\blacksquare$\\[2mm]} \setlength{\parindent}{0in} \def\E{\mathbb{E}} \def\Var{\text{Var}} \def\Cov{\text{Cov}} \def\Corr{\text{Corr}} \def\bW{\mathbf{W}} \newcommand{\bs}[1]{\bm{#1}} \newcommand{\Rho}{\mathrm{P}} \newcommand{\cif}{\text{ if }} \newcommand{\Ra}{\Longrightarrow} \newcommand{\ra}{\longrightarrow} \def\bphi{\boldsymbol \phi} \lstloadlanguages{R} \definecolor{keywordcolor}{rgb}{0,0.6,0.6} \definecolor{delimcolor}{rgb}{0.461,0.039,0.102} \definecolor{Rcommentcolor}{rgb}{0.101,0.043,0.432} \lstdefinestyle{Rsettings}{ basicstyle=\ttfamily, breaklines=true, showstringspaces=false, keywords={if, else, function, theFunction, tmp}, % Write as many keywords otherkeywords={}, commentstyle=\itshape\color{Rcommentcolor}, keywordstyle=\color{keywordcolor}, moredelim=[s][\color{delimcolor}]{"}{"}, } \lstset{basicstyle=\ttfamily, numbers=none, literate={~} {$\sim$}{2}} \def\ctitle{MATH 545: Introduction to Time Series Analysis} \def\ttitle{Winter 2025} \def\stitle{Conditional Heteroscedastic Models} \begin{document} <>= library(knitr) # global chunk options opts_chunk$set(cache=TRUE, autodep=TRUE) options(scipen=999) options(repos=c(CRAN="https://cloud.r-project.org/")) @ \maketitle In a \textit{conditional heteroscedastic model}, we have an observed data set with \textit{time-varying variance}. Typically, these models can be considered as special forms of state space models. We typically assume that \[ Y_t = \sigma_t Z_t \] where $\{Z_t\}$ is a Gaussian white noise process with variance 1; the volatility is therefore encapsulated in $\sigma_t$, and we adopt \begin{itemize} \item deterministic, observation-driven models, \item stochastic, parameter-driven models. \end{itemize} \bigskip \textbf{ARCH/GARCH MODELS:} The \text{ARCH} and \text{GARCH} models use \textit{deterministic} state space forms that allows the variance of the observed process to vary with time. The ARCH(1) model specifies $\sigma_t = \sqrt{X_t}$ so that \[ Y_t = \sqrt{X_t} Z_t \] and then that \[ X_{t} = \alpha_0 + \alpha_1 Y_{t-1}^2 \] whereas the GARCH(1,1) model specifies \[ X_{t} = \alpha_0 + \alpha_1 Y_{t-1}^2 + \beta_1 X_{t-1} \] where, in both cases, given $X_t = x_t$, $Y_t$ has mean zero and variance $x_t$. The most general form of the GARCH model is the $\textrm{GARCH}(p,q)$ model where \[ X_{t} = \alpha_0 + \sum_{j=1}^p \alpha_j Y_{t-j}^2 + \sum_{j=1}^q \beta_j X_{t-j}. \] \textbf{Moment properties:} For $\textrm{GARCH}(1,1)$, by standard properties, we have that $\E[Y_t] = 0$, and using iterated expectation \[ \Var[Y_t] = \E[Y_t^2] = \E[X_t Z_t^2] = \E[X_t] = \E[\alpha_0 + \alpha_1 Y_{t-1}^2 + \beta_1 X_{t-1}]. \] As $\Var[Y_t] = \E[Y_t^2] = \E[X_t] = \E[X_{t-1}] = \E[Y_{t-1}^2]$ for all $t$, we may conclude (by analogy with results for $\AR(1)$) that under stationarity \[ \Var[Y_t] = \alpha_0 + \alpha_1 \Var[Y_t] + \beta_1 \Var[Y_t] \] so that \[ \Var[Y_t] = \frac{1}{1-\alpha_1 - \beta_1} = \sigma_Y^2, \] say, and that a condition for stationarity is that $\alpha_1 + \beta_1 < 1$. It is evident also that for $h=\pm 1,\pm 2, \ldots$ \[ \E[Y_t Y_{t+h}] = \E[\sqrt{X_t} Z_t \sqrt{X_{t+h}} Z_{t+h}] = \E[\sqrt{X_t X_{t+h}}] \E[ Z_t Z_{t+h}] = 0 \] as $\{Z_t\}$ is white noise. However, \[ \Var[Y_t^2] = \E[Y_t^4] - \{\E[Y_t^2]\}^2 = \E[X_t^2 Z_t^4] - \sigma_Y^4 = 3 \E[X_t^2] - \sigma_Y^4 \] where \[ \E[X_t^2] = \Var[X_t] + \{E[X_t]\}^2 = \frac{\sigma_W^2}{1-\phi^2} + \mu_X^2 = \sigma_X^2 + \mu_X^2 \] say. Also, for $h=\pm 1, \pm 2, \ldots$, \[ \E[Y_t^2 Y_{t+h}^2] = \E[X_t Z_t^2 X_{t+h} Z_{t+h}^2] = \E[X_t X_{t+h}] \E[ Z_t^2 Z_{t+h}^2] = \E[X_t X_{t+h}] \neq 0 \] Then \begin{align*} \E[Y_t^2 Y_{t+h}^2] & = \E[(\alpha_0 + \alpha_1 Y_{t-1}^2 + \beta_1 X_{t-1})(\alpha_0 + \alpha_1 Y_{t+h-1}^2 + \beta_1 X_{t+h-1})]\\[6pt] & = \alpha_0^2 + \alpha_0 \alpha_1 \E[Y_{t+h-1}^2] + \alpha_0 \beta_1 \E[ X_{t+h-1}] \\[6pt] & \quad + \alpha_0 \alpha_1 \E[Y_{t-1}^2] + \alpha_1^2 \E[Y_{t-1}^2 Y_{t+h-1}^2] + \alpha_1 \beta_1 \E[X_{t+h-1} Y_{t-1}^2] \\[6pt] & \quad + \alpha_0 \beta_1 \E[X_{t-1}] + \alpha_1 \beta_1 \E[X_{t-1} Y_{t+h-1}^2] + \beta_1^2 \E[ X_{t-1} X_{t+h-1}] \end{align*} In this expression, \begin{align*} \E[X_t X_{t+h}] & = \gamma_{Y^2}(h)\\ \E[Y_{t+h-1}^2] & = \sigma_Y^2 \\ \E[ X_{t+h-1}] & = \alpha_0 + \alpha_1 \E[Y_{t+h-2}^2] = \alpha_0 + \alpha_1 \sigma_Y^2\\ \E[Y_{t-1}^2 Y_{t+h-1}^2] & = \gamma_{Y^2}(h)\\ \E[X_{t+h-1} Y_{t-1}^2] & = \E[(\alpha_0 + \alpha_1 Y_{t+h-2}^2 ) Y_{t-1}^2 ] = \alpha_0 \sigma_Y^2 + \alpha_1 \gamma_{Y^2}(h-1)\\ \E[X_{t-1} Y_{t+h-1}^2] & = \E[(\alpha_0 + \alpha_1 Y_{t-2}^2 ) Y_{t+h-1}^2 ] =\alpha_0 \sigma_Y^2 + \alpha_1 \gamma_{Y^2}(h+1)\\ \E[ X_{t-1} X_{t+h-1}] & = \gamma_{Y^2}(h) \end{align*} so \begin{align*} \E[Y_t^2 Y_{t+h}^2] & = \alpha_0^2 + \alpha_0 \alpha_1 \sigma_Y^2 + \alpha_0 \beta_1 (\alpha_0 + \alpha_1 \sigma_Y^2)\\[6pt] & \quad + \alpha_0 \alpha_1 \sigma_Y^2 + \alpha_1^2 \gamma_{Y^2}(h) + \alpha_1 \beta_1 (\alpha_0 \sigma_Y^2 + \alpha_1 \gamma_{Y^2}(h-1)) \\[6pt] & \quad + \alpha_0 \beta_1 (\alpha_0 + \alpha_1 \sigma_Y^2) + \alpha_1 \beta_1 (\alpha_0 \sigma_Y^2 + \alpha_1 \gamma_{Y^2}(h+1)) )+ \beta_1^2 \gamma_{Y^2}(h) \end{align*} This ultimately gives a linear recursion of the form \[ \gamma_{Y^2}(h+1) = A + B \gamma_{Y^2}(h) + C \gamma_{Y^2}(h-1) \] with initial conditions obtained by examining $h=0,1$. \bigskip \textbf{STOCHASTIC VOLATILITY MODELS:} The \textit{stochastic volatility} model uses a stochastic state space form that allows the variance of the observed process to vary stochastically with time. The simplest form specifies that \[ Y_t = \exp\{X_t/2\} Z_t \qquad \{Z_t\} \sim \WN(0,\sigma_Z^2) \] and then that \[ X_{t+1} = \phi X_t + W_{t} \qquad \{W_t\} \sim \WN(0,\sigma_W^2) \] that is, given $X_t = x_t$, $Y_t$ has mean zero and variance $\exp\{x_t\}$, and $\{X_t\} \sim \AR(1)$. We can set $\sigma_Z^2=1$, and can reparameterize the model as follows: \[ Y_t = \exp\{X_t/2\} Z_t \] with $\{Z_t\} \sim \WN(0,1)$, and \[ X_{t+1} = \phi_0 + \phi X_t + W_{t} \] or equivalently \[ X_{t+1} = \mu_X + \phi(X_t-\mu_X) + W_{t} \] where $\phi_0 = (1-\phi) \mu_X $, and $\{W_t\} \sim \WN(0,\sigma_W^2)$. For this model it is evident that $\E[Y_t] = 0$ for all $t$, and \[ \Var[Y_t] = \E[Y_t^2] = \E[\exp\{X_t\} Z_t^2] = \E[\exp\{X_t\}], \] whereas under stationarity, if $\{Z_t\}$ and $\{W_t\}$ are uncorrelated \[ \gamma_Y(h) = \E[Y_t Y_{t+h}] = \E[\exp\{X_t/2\} Z_t \exp\{X_{t+h}/2\} Z_{t+h}] = \E[\exp\{(X_t+X_{t+h})/2\}] \E[Z_t Z_{t+h}] = 0. \] Now \[ \Var[Y_t^2] = \E[Y_t^4] - \{\E[Y_t^2]\}^2 = \E[\exp\{2 X_t\} Z_t^4] - \{\E[\exp\{X_t\}]\}^2 = \E[\exp\{2 X_t\}] \E[Z_t^4] - \{\E[\exp\{X_t\}]\}^2. \] Under the assumption of Normality of the $\{Z_t\}$ process, we have $\E[Z_t^4] = 3$ by properties of the standard Normal. Similarly \[ \E[Y_t^2 Y_{t+h}^2] = \E[\exp\{X_t\} Z_t^2 \exp\{X_{t+h}\} Z_{t+h}^2] = \E[\exp\{(X_t+X_{t+h})\}] \E[Z_t^2 Z_{t+h}^2] = \E[\exp\{(X_t+X_{t+h})\}] \] as $\E[Z_t^2] = 1$. So for $h \geq 1$, \[ \gamma_{Y^2}(h) = \E[Y_t^2 Y_{t+h}^2] - \{\E[Y_t^2]\}^2 = \E[\exp\{(X_t+X_{t+h})\}] - \{\E[\exp\{X_t\}]\}^2 \] \bigskip To study properties of this model we use properties of the log-normal distribution that allow us to compute \[ \E[\exp\{X_t\}] \qquad \E[\exp\{2 X_t\}] \qquad \E[\exp\{(X_t+X_{t+h})\}]. \] Specifically \[ \E \left[ \exp \left( X_{t}\right) ^{j}\right] =\exp \left\{ j\mu _{X}+\frac{1}{2} j^{2}\sigma_{X}^{2}\right\}, \] so that, if $r = 2 s$ is even and $\{X_{t}\}$ is stationary, all the even moments of $Y_{t}$ exist and are given by the formula \begin{equation*} \E \left[ Y_{t}^{r}\right] =\E \left[ \left( \exp (X_{t}) \right) ^{s}\right] \E \left[ Z^{r}\right] =\exp \left\{ s \mu_{X} +\dfrac{s^{2}}{2}\sigma_{X}^{2}\right\} \frac{(2s)!}{2^{s}\ s !}. \label{Nrthmom} \end{equation*} whereas all the odd moments are zero. Thus \begin{align*} \Var\left[ Y_{t} \right] & =\E \left[ Y_{t}^{2}\right] = \E[\exp\{X_t\}] =\exp \left\{ \mu_{X}+ \frac{1}{2}\sigma_{X}^{2}\right\} \\[6pt] \Var\left[ Y_{t}^2 \right] & = 3 \E[\exp\{2 X_t\}] - \{\E[\exp\{X_t\}]\}^2 = 3 \exp\left\{2 \mu_X + 2 \sigma_X^2 \right\} - \exp\left\{2 \mu_X + \sigma_X^2 \right\}. \end{align*} Now, for this model, we know that $X_t + X_{t+h}$ is a Normal random variable with mean $2 \mu_X$. We have by standard arguments that \begin{align*} \Cov[X_t,X_{t+h}] & = \Cov[X_t,\mu_X + \phi( X_{t+h-1} - \mu)] \\[6pt] & = \Cov[X_t,\phi X_{t+h-1}] \\[6pt] & = \Cov[X_t, \phi^h X_t] & \text{recursion}\\[6pt] & = \phi^h \Var[X_t] \\[6pt] & = \phi^h \sigma_X^2. \end{align*} Therefore we have that \[ \Var[X_t + X_{t+h}] = \Var[X_t] + \Var[X_{t+h}] + 2 \Cov [X_t,X_{t+h}] = 2\sigma_X^2 + 2 \phi^h \sigma_X^2. \] so that \[ X_t + X_{t+h} \sim Normal(2 \mu_X, 2 \sigma_X^2 (1+\phi^{h})). \] and hence for $h \geq 1$ \[ \gamma_{Y^2}(h) = \E[\exp\{(X_t+X_{t+h})\}] - \{\E[\exp\{X_t\}]\}^2 = \exp\left \{ 2 \mu_X + \sigma_X^2 (1+\phi^{h}) \right\} - \exp\left \{ 2 \mu_X + \sigma_X^2 \right\} \] which shows that $Y_t^2$ is an autocorrelated process. We have that the ACF for $h = \pm 1, \pm 2,\ldots$ is \begin{align*} \rho_{Y^2}(h) & = \frac{\exp\left \{ 2 \mu_X + \sigma_X^2 (1+\phi^{h}) \right\} - \exp\left \{ 2 \mu_X + \sigma_X^2 \right\}}{3 \exp\left\{2 \mu_X + 2 \sigma_X^2 \right\} - \exp\left\{2 \mu_X + \sigma_X^2 \right\}} \\[6pt] & = \frac{\exp\left \{\sigma_X^2 \phi^{h}) \right\} - 1}{3 \exp\left\{ \sigma_X^2 \right\} - 1} \end{align*} \pagebreak \textbf{Simulation study:} We begin with the GARCH models, and simulate from the ARCH(1) model with $\alpha_0 = 0.1$ and $\alpha_1 = 0.95$. <>= set.seed(3245) n<-1000;al0<-0.01;al1<-0.95 X<-Y<-rep(0,n);X[1]<-al0;Y[1]<-rnorm(1)*sqrt(X[1]) #ARCH(1) MODEL for(j in 2:n){ X[j]<-al0+al1*Y[j-1]^2 Y[j]<-rnorm(1)*sqrt(X[j]) } @ GARCH models can be fitted using the \texttt{fGarch} library. <>= library(fGarch) mydata<-data.frame(Y) #Fitting the ARCH(1) model fit.10<-garchFit(Y~garch(1,0),data=mydata,trace=F) #ARCH(1) show(fit.10) #Fitting the GARCH(1,1) model fit.11<-garchFit(Y~garch(1,1),data=mydata,trace=F) #GARCH(1,1) show(fit.11) c(fit.10@fit$ics[1],fit.11@fit$ics[1]) #AIC comparison fit.10@fit$coef fit.11@fit$coef @ Thus the maximum likelihood estimated model is \begin{align*} \text{ARCH}(1) \qquad X_t & = \Sexpr{fit.10@fit$coef[2]} + \Sexpr{fit.10@fit$coef[3]} Y_{t-1}^2\\ \text{GARCH}(1,1) \qquad X_{t} & = \Sexpr{fit.11@fit$coef[2]} + \Sexpr{fit.11@fit$coef[3]} Y_{t-1}^2 + \Sexpr{fit.11@fit$coef[4]} X_{t-1} \end{align*} The ARCH(1) model is preferred under the AIC comparison. The parameter referred to as \texttt{omega} is $\alpha_0$, and \texttt{mu} is the mean of the $\{Y_t\}$ series. If we inspect the residual plot resulting from the ARCH(1) fit, the variance of the residuals seem to be adequate, apart from in one part of the plot. <>= plot(residuals(fit.10),type='l',xlab='t',ylab='Residuals') @ \bigskip \textbf{Stochastic volatility simulation:} For the SV simulation, in the original parameterization we have \[ Y_t = \exp\{X_t/2\} Z_t \] with $\sigma_Z^2 = 0.01$, and \[ X_t = \phi X_{t-1} + W_t \] with $\sigma_W^2 = 0.05$. If we re-scale $Z_t$ to be unit variance, then this the same model as obtained by re-writing the state equation as \[ X_t = \mu_X + \phi(X_{t-1} - \mu_X) + W_t \] with $\mu_X = 2 \log \sigma_Z^2$. <>= set.seed(24) n<-1000;phi<-0.98;sigZ<-0.01;sigW<-0.05 X<-arima.sim(n,model=list(ar=phi),rand.gen = function(n) sqrt(sigW) * rnorm(n)) Z<-rnorm(n)*sigZ Y0<-exp(X/2)*Z set.seed(24) mu<-2*log(sigZ) #Reparameterization X<-arima.sim(n,model=list(ar=phi),rand.gen = function(n) sqrt(sigW) * rnorm(n)) X<-X+mu Z<-rnorm(n) Y<-exp(X/2)*Z par(mar=c(4,3,2,0));plot(100*Y0,type='l',ylim=range(-3.5,3.5)); points(100*Y0,pch=19,cex=0.3) plot(100*Y,type='l',ylim=range(-3.5,3.5)); points(100*Y,pch=19,cex=0.3) @ We can compare the theoretical ACF with the sample ACF <>= muX<-mu sigX<-sigW/(1-phi^2) gamma0<-3*exp(2*muX+2*sigX)-exp(2*muX+sigX) hvec<-c(1:30) gammah<-exp(2*muX+sigX*(1+phi^hvec))-exp(2*muX+sigX) gammah/gamma0 #ACF acf.vec<-c(gamma0,gammah)/gamma0 acf(Y^2,lag.max=30,,main=expression(paste('ACF of ',Y[t]^2))) points(c(0,hvec),acf.vec,pch=19,cex=0.5,col='red') legend(15,1.0,c('True ACF'),pch=19,col='red') @ \pagebreak Inference for the SV model can be carried out using Bayesian methods and Markov chain Monte Carlo (MCMC) and the \texttt{stochvol} library. \begin{itemize} \item \texttt{svsim} simulates a SV model \item \texttt{svsample} carries out MCMC analysis \end{itemize} <>= library(stochvol) y.sim <- svsim(1000, mu = 2*log(sigZ), phi = phi, sigma = sigW) plot(y.sim) @ \medskip The MCMC algorithm draws 5000 samples from the \textit{posterior distribution} for the three parameters of the model $\mu,\phi$ and $\sigma_W$. <>= mcmc.draws <- svsample(y.sim$y, draws = 5000, burnin = 1000, priormu = c(0, 10), priorphi = c(20, 1.2), priorsigma = 100) para<-data.matrix(mcmc.draws$para) par(mfrow=c(2,2),mar=c(4,4,2,0)) hist(para[,1],col='gray',nclass=30,xlab='', main=expression(paste('Posterior for ',mu))) box();abline(v=mu,col='red',lwd=2) #Mark the true values with red lines. hist(para[,2],col='gray',nclass=30,xlab='', main=expression(paste('Posterior for ',phi))) box();abline(v=phi,col='red',lwd=2) hist(para[,3],col='gray',nclass=30,xlab='', main=expression(paste('Posterior for ',sigma[W]^2))) box();abline(v=sigW,col='red',lwd=2) @ Parameter estimates can be obtained from the posterior distribution by summarizing the posterior samples: <>= apply(para,2,summary) @ We can also forecast using the model, with commands <>= fore <- predict(mcmc.draws, 20) plot(mcmc.draws, forecast = fore) @ \pagebreak \textbf{Real Data Analysis:} For an analysis of a real data set, we consider the S \& P 500 index over the period 2007-2017. <>= SP.2017<-read.csv('GSPC.csv',header=TRUE) X<-SP.2017$Close Tv<-as.Date(as.character(SP.2017$Date),'%d/%m/%Y') par(mar=c(4,4,2,0)) plot(Tv,X,type='l',main='S & P 500 index',xlab='Date') @ If we let $X_t$ denote the index value, we then transform to the percentage log returns series \[ Y_t = 100 \log (X_{t+1}/X_t) = 100(\log X_{t+1} - \log X_{t}) \] <>= Y<-100*diff(log(X)) par(mar=c(4,4,2,0)) plot(Tv[-1],Y,type='l',main='S & P 500 log returns',xlab='Date') hist(Y,nclass=50,col='gray',main='Histogram of log returns');box() @ The histogram indicates a \textit{leptokurtic} (or heavy-tailed) marginal distribution for the log returns. <>= par(mar=c(4,4,4,0)) acf(Y,lag.max=100,main=expression(paste('ACF of ',Y[t]))) acf(Y^2,lag.max=100,main=expression(paste('ACF of ',Y[t]^2))) @ \medskip The log returns series is uncorrelated, but the squared series has significant autocorrelation. \pagebreak The ARCH(1) and the GARCH(1,1) models can be used. <>= mydata<-data.frame(Y) fit.10<-garchFit(Y~garch(1,0),data=mydata,trace=F) #ARCH(1) show(fit.10) fit.11<-garchFit(Y~garch(1,1),data=mydata,trace=F) #GARCH(1,1) show(fit.11) c(fit.10@fit$ics[1],fit.11@fit$ics[1]) #AIC comparison @ Here the GARCH(1,1) model is preferred, and we have \[ X_{t} = \Sexpr{fit.11@fit$coef[2]} + \Sexpr{fit.11@fit$coef[3]} Y_{t-1}^2 + \Sexpr{fit.11@fit$coef[4]} X_{t-1} \] as the maximum likelihood fitted model. Here we have \[ \widehat \alpha_1 + \widehat \beta_1 = \Sexpr{fit.11@fit$coef[3]} + \Sexpr{fit.11@fit$coef[4]} = \Sexpr{fit.11@fit$coef[3]+fit.11@fit$coef[4]} < 1 \] \bigskip We can carry out an analysis using the stochastic volatility model. Here we use the MCMC functionality in the \texttt{stochvol} package. <>= mcmc.draws <- svsample(Y, draws = 5000, burnin = 1000, priormu = c(0, 10), priorphi = c(20, 1.2), priorsigma = 100) para<-data.matrix(mcmc.draws$para) par(mfrow=c(2,2),mar=c(4,4,2,0)) hist(para[,1],col='gray',nclass=30,xlab='', main=expression(paste('Posterior for ',mu))) box() hist(para[,2],col='gray',nclass=30,xlab='', main=expression(paste('Posterior for ',phi))) box() hist(para[,3],col='gray',nclass=30,xlab='', main=expression(paste('Posterior for ',sigma[W]^2))) box() @ <>= fore <- predict(mcmc.draws, 20) plot(mcmc.draws, forecast = fore) @ \end{document}