\documentclass[notitlepage]{article} \usepackage{../../Math545} \usepackage{listings} \usepackage{numprint} \usepackage{enumerate} \usepackage{bbm} \usepackage{bm} \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\phivec{\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{Model-based Estimation for ARMA models} \begin{document} <>= library(knitr) # global chunk options opts_chunk$set(cache=TRUE, autodep=TRUE) options(scipen=999) @ \maketitle For a zero mean $\ARMA(p,q)$ model, we can make most efficient estimation of the autocovariance and most accurate forecasting using \textit{model-base}d estimation provided we can correctly specify the model concerned. This is because the acvf is entirely specified by the $p+q+1$ parameters $\phi_1,\ldots,\phi_p, \theta_1,\ldots,\theta_q,\sigma_Z^2$, and typically $p+q+1$ is much smaller than $n$. \bigskip \textbf{Yule-Walker Estimation:} The Yule-Walker method can be applied to zero mean, stationary $\AR(p)$ processes; \[ X_t = \sum_{j=1}^p \phi_j X_{t-j} + Z_t. \] We multiply both sides by $X_{t-h}$ for $h=1,\ldots,p$, and take expectations to get \[ \E[X_t X_{t-h}] = \sum_{j=1}^p \phi_j \E[X_{t-j} X_{t-h}] \qquad \Longrightarrow \qquad \gamma_X(h) = \sum_{j=1}^p \phi_j \gamma_X(h-j) \] for $h=1,\ldots,p$. Taking these $p$ equations together, we obtain $\Gammat_p \phivec = \gammavec_p$, that is \[ \begin{bmatrix} \gamma_X(0) & \gamma_X(1) & \gamma_X(2) & \cdots & \gamma_X(p-1) \\[6pt] \gamma_X(1) & \gamma_X(0) & \gamma_X(1) & \cdots & \gamma_X(p-2) \\[6pt] \gamma_X(2) & \gamma_X(1) & \gamma_X(0) & \cdots & \gamma_X(p-3) \\[6pt] \vdots & \vdots & \vdots & \ddots & \vdots\\[6pt] \gamma_X(p-1) & \gamma_X(p-2) & \gamma_X(p-3) & \cdots & \gamma_X(0) \end{bmatrix} \begin{bmatrix}\phi_1 \\[6pt] \phi_2 \\[6pt] \phi_3 \\[6pt] \vdots \\[6pt] \phi_p \end{bmatrix}= \begin{bmatrix}\gamma_X(1) \\[6pt] \gamma_X(2) \\[6pt] \gamma_X(3) \\[6pt] \vdots \\[6pt] \gamma_X(p)\end{bmatrix}. \] We may now substitute in the empirical estimates for the acv values in the matrix and in the vector on the right-hand side, and solve the resulting equation to obtain the estimates $\widehat {\phivec} = (\widehat \phi_1, \ldots, \widehat \phi_p)^\top$. \medskip \textbf{EXAMPLE:} Consider the $\AR(4)$ model given by $\Phi(B) X_t = Z_t$, where polynomial $\Phi(\cdot)$ is given by \[ \Phi(z) = (1-0.5 z) (1+0.4z) (1-0.2z)^2 \] so that the reciprocal roots are $0.5,-0.4,0.2$ and $0.2$. This means that the model is stationary and causal, and the coefficients are computed as follows <>= alpha.vals<-c(0.5,-0.4,0.2,0.2) ar.coefs<-rep(1,length(alpha.vals)+1) for(j in 1:length(alpha.vals)){ ar.coefs.old<-ar.coefs if(j > 1){ for(k in 2:j){ ar.coefs[k]<-ar.coefs.old[k]-alpha.vals[j]*ar.coefs.old[k-1] } } ar.coefs[j+1]<--alpha.vals[j]*ar.coefs.old[j] } (phi.vec<--ar.coefs[-1]) #These are the AR coefficients @ that is $\phi_1 = \Sexpr{phi.vec[1]}, \phi_2 = \Sexpr{phi.vec[2]} ,\phi_3 = \Sexpr{phi.vec[3]},\phi_4 = \Sexpr{phi.vec[4]}$. Numerically we have <>= p<-4 gammaX<-ARMAacf(ar=phi.vec,lag.max=p);gammaX gammaX.mat<-toeplitz(gammaX[-(p+1)]) solve(gammaX.mat,gammaX[-1]) #Check that the true phi values are a solution @ Now we repeat the calculation with a simulated data set, replacing $\Gammat_p$ by $\widehat \Gammat_p$. <>= n<-500 X<-arima.sim(n,model=list(ar=phi.vec)) gammaX.hat<-acf(X,plot=F,lag.max=p)$acf gammaX.mat.hat<-toeplitz(gammaX.hat[-(p+1)]) solve(gammaX.mat.hat,gammaX.hat[-1]) #Estimates @ These estimates can be obtained using the \texttt{ar} function: <>= ar(X,order.max=4,aic=FALSE,method='yule-walker') @ \bigskip \textbf{Ordinary Least Squares:} We can also obtain parameter estimates for an $\AR$ model by ordinary least squares (OLS) by minimizing \[ \sum_{t=p+1}^n \left(x_t - \sum_{j=1}^p \phi_j x_{t-j} \right)^2 \] which is a straightforward problem that can be solved analytically using calculus. In \texttt{R} to solve the OLS problem we may use the \texttt{lm} function. In the example above, if we want to assume that the process is zero mean, then we proceed as follows: <>= Y<-X[(p+1):n] X1<-X[p:(n-1)] X2<-X[(p-1):(n-2)] X3<-X[(p-2):(n-3)] X4<-X[(p-3):(n-4)] summary(lm(Y~-1+X1+X2+X3+X4)) @ The \texttt{ar} function can also be used to obtain the estimates. <>= ar(X,order.max=4,demean=FALSE,aic=FALSE,method='ols') @ We can also proceed by estimating the mean and subtracting it prior to the estimation: <>= mu.hat<-mean(X) mu.hat #Estimated mean Y<-X[(p+1):n]-mu.hat X1<-X[p:(n-1)]-mu.hat X2<-X[(p-1):(n-2)]-mu.hat X3<-X[(p-2):(n-3)]-mu.hat X4<-X[(p-3):(n-4)]-mu.hat summary(lm(Y~-1+X1+X2+X3+X4)) ar(X-mu.hat,order.max=4,demean=FALSE,aic=FALSE,method='ols') @ \bigskip \textbf{Burg's method:} Burg's method for an $\AR$ process uses the partial autocorrelations to lag $p$. Using the \texttt{ar} function: <>= ar(X,order.max=4,demean=FALSE,aic=FALSE,method='burg') @ \textbf{Maximum likelihood:} A maximum likelihood approach requires a further assumption that specifies the conditional probability distribution of each of the $X_t$. Usually it is assumed that $Z_t \sim Normal(0,\sigma_Z^2)$, so that for $t > p$ \[ X_{t}|X_1,\ldots,X_{t-1} \sim Normal(\phi_1 X_{t-1} + \phi_2 X_{t-2} + \cdots + \phi_p X_{t-p}, \sigma_Z^2) \] from which we can construct the likelihood function in the usual way. For $t = 1,\ldots,p$, we may deduce the distribution of $X_t$ using the stationarity of the process, and deduce that \[ X_{1:n} \sim Normal_n (\mathbf{0},\Gammat_n) \] Here $\Gammat(n)=\Gammat_n(\phivec)$, so that the full likelihood for data $\bx_{1:n}$ is given by \begin{equation} \label{eq:like} \mathcal{L}_n(\phivec) = f(x_1,\ldots,x_p;\phivec) \prod_{t=p+1}^n f(x_t|x_1,\ldots,x_{t-1};\phivec). \end{equation} Up to a constant the corresponding log-likelihood is \[ \ell_n(\phivec) = - \frac{1}{2} \log | \Gammat_n(\phivec)| - \frac{1}{2} \bx_{1:n}^\top \left\{\Gammat_n(\phivec)\right\}^{-1} \bx_{1:n}. \] To maximize this function, we need to compute the matrix inverse and the determinant. This can be achieved efficiently using the Levinson-Durbin or Innovations algorithms, or using eigenvalue/singular value decompositions where we write \[ \Gammat_X(n;\phivec) = \bU \bD \bU^\top \qquad \therefore \qquad \{\Gammat_X(n;\phivec)\}^{-1} = \bU \bD^{-1} \bU^\top \] with $\bD$ a diagonal matrix with elements $d_1,\ldots,d_n$, and \[ | \Gammat_n(\phivec)| \equiv \text{det} \ \Gammat_n(\phivec) = \prod_{t=1}^n d_t \] <<,comment='+', fig.width=6.25, fig.height=4.0, fig.align='center'>>= ar(X,order.max=4,demean=FALSE,aic=FALSE,method='mle') @ The maximum likelihood method can be used for $\ARMA(p,q)$ processes. \medskip \textbf{EXAMPLE:} Consider the $\ARMA(2,3)$ model \[ X_t = 0.1 X_{t-1} - 0.56 X_{t-2} + Z_t + 0.2 Z_{t-1} -0.53 Z_{t-2} + 0.09 Z_{t-3}. \] In this model, we have \[ \Phi(z) = 1 - 0.1 z - 0.56 z^2 = (1-0.8z)(1+0.7z) \] so the model is stationary and causal, and \[ \Theta(z) = (1-0.2z) (1-0.5 z) (1+0.9z) \] so the model is invertible. The $\Theta(z)$ coefficients are computed as <>= beta.vals<-c(0.2,0.5,-0.9) ma.coefs<-rep(1,length(beta.vals)+1) for(j in 1:length(beta.vals)){ ma.coefs.old<-ma.coefs if(j > 1){ for(k in 2:j){ ma.coefs[k]<-ma.coefs.old[k]-beta.vals[j]*ma.coefs.old[k-1] } } ma.coefs[j+1]<--beta.vals[j]*ma.coefs.old[j] } (theta.vec<-ma.coefs[-1]) @ The acf of this model shows slow decay <>= phi.vec<-c(0.1,0.56);true.pars<-c(phi.vec,theta.vec) acfs<-ARMAacf(ar=phi.vec,ma=theta.vec,lag.max=50) acfs[1:10] par(mar=c(4,4,1,0)) plot(0:50,acfs,pch=19,cex=0.5,xlab='Lag',ylab='True acf') for(i in 0:50){lines(c(i,i),c(0,acfs[1+i]))} @ \medskip We start with a very large sample of $n=2000$: <>= n<-2000 set.seed(333) X<-arima.sim(n,model=list(ar=phi.vec,ma=theta.vec)) (ml1<-arima(X,order=c(2,0,3),method='ML',include.mean = FALSE)) @ The estimation procedure recovers the true values of the parameters quite accurately. For a smaller sample size, the estimation is more inaccurate. <>= n<-200 set.seed(333) X<-arima.sim(n,model=list(ar=phi.vec,ma=theta.vec)) (ml2<-arima(X,order=c(2,0,3),method='ML',include.mean = FALSE)) @ We can run a simulation study to assess the accuracy of inference for the different sample size, and run 1000 replications. The following code runs the simulation and stores the parameter estimates. <>= n1<-2000;n2<-200;nreps<-1000 ests.mat1<-ests.mat2<-matrix(0,nrow=nreps,ncol=5) oldw <- getOption("warn") options(warn=-1) for(irep in 1:nreps){ X<-arima.sim(n1,model=list(ar=phi.vec,ma=theta.vec)) ml1<-arima(X,order=c(2,0,3),method='ML',include.mean = FALSE) ests.mat1[irep,]<-coef(ml1) X<-arima.sim(n2,model=list(ar=phi.vec,ma=theta.vec)) ml2<-arima(X,order=c(2,0,3),method='ML',include.mean = FALSE) ests.mat2[irep,]<-coef(ml2) } options(warn=oldw) par(mfrow=c(2,3),mar=c(4,4,2,0),oma=c(0,0,1,0)) tvec<-c(expression(hat(phi)[1]),expression(hat(phi)[2]), expression(hat(theta)[1]),expression(hat(theta)[2]),expression(hat(theta)[3])) for(i in 1:5){ hist(ests.mat1[,i],col='gray',main=tvec[i],xlab=' ');abline(v=true.pars[i],col='red') box() } title(substitute(paste('n =',nval),list(nval = n1)), outer=TRUE) par(mfrow=c(2,3),mar=c(4,4,2,0),oma=c(0,0,1,0)) for(i in 1:5){ hist(ests.mat2[,i],col='gray',main=tvec[i],xlab=' ');abline(v=true.pars[i],col='red') box() } title(substitute(paste('n =',nval),list(nval = n2)), outer=TRUE) @ The red lines mark the positions of the true values: \bigskip \textbf{Note:} In all of the inference methods, it is the information in the autocovariances that determines quality of the parameter estimation. As different $\ARMA$ models can have almost identical acvs, it may be difficult in general to identify particular parameters. \medskip \textbf{EXAMPLE:} Consider the $\AR(2)$ model \[ X_t = X_{t-1} - 0.09 X_{t-2} + Z_t \] and the $\ARMA(3,2)$ model \[ X_t = 1.1 X_{t-1} - 0.19 X_{t-2} + 0.009 X_{t-3} + Z_t - 0.25 Z_{t-1} + 0.2 Z_{t-2} \] <>= phi.vec<-c(1,-0.09) acf.1<-ARMAacf(ar=phi.vec,lag.max=20) phi.vec<-c(1.1,-0.19,0.009);theta.vec<-c(-0.25,0.2) acf.2<-ARMAacf(ar=phi.vec,ma=theta.vec,lag.max=20) par(mar=c(4,4,1,0)) plot(0:20,acf.1,pch=19,cex=0.5,xlab='Lag',ylab='True acf',col='red',ylim=range(0,1)) points(0:20,acf.2,pch=19,cex=0.5,col='blue') legend(10,1.0,c("Model 1","Model 2"),pch=19,col=c('red','blue')) @ Clearly the acf values are almost identical. \end{document}