\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} \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{Autocovariance Estimation for Time Series Data} \begin{document} <>= library(knitr) # global chunk options opts_chunk$set(cache=TRUE, autodep=TRUE) options(scipen=999) @ \maketitle When presented with a real time series data set $\bx_{1:n} = (x_1,\ldots,x_n)$ it is of interest to attempt to learn which model might be appropriate for modelling the variation. We focus on stationary models -- this assumes that appropriate pre-processing in the form of \begin{itemize} \item de-meaning (eg trend removal) \item removal of seasonal (periodic) components \item differencing or seasonal differencing \end{itemize} has been carried out prior to the analysis, so that we may assume that the data can be represented by a stationary, zero-mean process. \bigskip \textbf{Autocovariance Estimation:} A common first task is to identify the autocovariance (and autocorrelation) structure. As we have seen previously, the recommended moment-based (or model-free) estimate of the autocovariance at lag $h$ is \[ \widehat \gamma_X(h) = \frac{1}{n} \sum_{t=1}^{n-h} (x_{t+h} - \overline{x}) (x_{t} - \overline{x}) \qquad h = 0,1,\ldots,h_{\max} \] where $h_{\max}$ is chosen such that the corresponding estimator has reasonable properties; in practice, we cannot choose $h_{\max}$ too large as this reduces the number of terms in the sum by too large an amount, which renders the estimator highly variable. By standard statistical theory, we have that as $n \longrightarrow \infty$ that for the estimator \[ \sqrt{n}(\widehat \gamma_X(h) - \gamma_X(h)) \CiD Normal(0,v_X(h)) \qquad h \in \mathbb{Z} \] for some \textit{asymptotic variance} $v_X(h)$. However, note that for a finite data set the corresponding large sample approximation \[ \widehat \gamma_X(h) \approxsim Normal(\gamma_X(h),\widehat v_X(h)/n) \] may not hold if $h$ is too large. \bigskip \textbf{EXAMPLE:} To illustrate the challenges, consider the $\ARMA(2,3)$ model \[ X_t = 0.5 X_{t-1} + 0.2 X_{t-2} + Z_t + Z_{t-1} + 2 Z_{t-2} - Z_{t-3} \] <>= phi.vec<-c(0.5,0.2);theta.vec<-c(1,2,-1);sigmaZ<-1 true.var<-sigmaZ^2*(1+sum(ARMAtoMA(ar=phi.vec,ma=theta.vec,lag.max=2000)^2)) true.var true.acvs<-true.var*ARMAacf(ar=phi.vec,ma=theta.vec,lag.max=20) true.acvs[1:11] true.acf<-true.acvs/true.var true.acf[1:11] @ An example realization from this process (assuming a Gaussian white noise process with variance $\sigma_Z^2=1$) can be seen below <>= set.seed(237);n<-200 X<-arima.sim(n,model=list(ar=phi.vec,ma=theta.vec)) par(mar=c(4,4,1,2));plot(X,type='l',xlab='Time');points(X,pch=19,cex=0.6) @ The \texttt{acf} function plots the sample acf: <>= par(mar=c(4,4,1,1)) acf(X,main='ACF for sample data',ylim=range(-1,1),lag.max=20) points(0:20,true.acf,pch=4,col='red') legend(14,1,c('True acf'),pch=4,col='red') @ The dotted blue horizontal lines are at $\pm 1.96/\sqrt{n}$; they represent confidence bounds for the correlation assuming a null model that the process is a white noise process. \bigskip We can examine the statistical performance of the moment-based estimator in a simulation (or Monte Carlo) study. We generate 5000 independent instances of the series, estimate the correlation for each, and then summarize the results in a boxplot of the `standardized' quantity \[ \sqrt{n}(\widehat \gamma_X(h) - \gamma_X(h)) \] We first do this for a very large sample of $n=2000$: <>= set.seed(237);n<-2000;nits<-5000;lmax<-20 gamma.hat<-rho.hat<-matrix(0,nrow=nits,ncol=lmax+1) for(iter in 1:nits){ X<-arima.sim(n,model=list(ar=phi.vec,ma=theta.vec)) gamma.hat[iter,]<-acf(X,lag.max=lmax,plot=F,type='covariance')$acf rho.hat[iter,]<-gamma.hat[iter,]/gamma.hat[iter,1] } #Standardize the estimates of the acv and the acf std.gamma.hat<-t(apply(gamma.hat,1,function(x){return(sqrt(n)*(x-true.acvs))})) std.rho.hat<-t(apply(rho.hat,1,function(x){return(sqrt(n)*(x-true.acf))})) par(mar=c(4,4,0,2)) boxplot(std.gamma.hat,col='gray',pch=19,cex=0.7);abline(h=0,lty=2,col='red') apply(std.gamma.hat,2,var) #Compute the sampling variance of the estimator boxplot(std.rho.hat,col='gray',pch=19,cex=0.7); abline(h=0,lty=2,col='red') apply(std.rho.hat,2,mean) #Compute the expected value of the estimator apply(std.rho.hat,2,var) #Compute the sampling variance of the estimator @ The Monte Carlo study suggests that the standardized bias defined as \[ \sqrt{n}\E[(\widehat \gamma_X(h) - \gamma_X(h))] \] is relatively low. It can be shown that the distribution of the vector of estimators of the first $h \ll n$ ACF values, $\boldsymbol{\rho}_X(1:h) = (\rho_X(1),\ldots,\rho_X(h))^\top$ is approximately multivariate Normal, that is, as $n \longrightarrow \infty$, we have that \[ \sqrt{n} (\widehat{\mathbf{\rho}}_X(1:h) - {\mathbf{\rho}}_X(1:h)) \CiD Normal_h(\mathbf{0},\bW) \] where $\bW$ is an $h \times h$ matrix with $(i,j)$th element \[ w_{ij} = \sum_{k=1}^\infty \{\rho_X(k+i) + \rho_X(k-i)-2\rho_X(i)\rho_X(k)\} \{\rho_X(k + j) + \rho_X(k-j)-2\rho_X(j)\rho_X(k)\}; \] this is known as \textbf{Bartlett's Formula}. Thus the diagonal elements are \[ w_{ii} = \sum_{k=1}^\infty \{\rho_X(k+i) + \rho_X(k-i)-2\rho_X(i)\rho_X(k)\}^2 \qquad i=1,\ldots,h. \] The matrix $\bW$ may be estimated by using the sample ACF values. <>= hmax<-20 W<-matrix(0,hmax,hmax) lmax<-2000 true.acf<-ARMAacf(ar=phi.vec,ma=theta.vec,lag.max=lmax+hmax) for(i in 1:hmax){ i1<-c(1:lmax)+i i2<-c(1:lmax)-i ti<-true.acf[i1+1]+true.acf[abs(i2)+1]-2*true.acf[i+1]*true.acf[1:lmax+1] for(j in 1:hmax){ j1<-c(1:lmax)+j j2<-c(1:lmax)-j tj<-true.acf[j1+1]+true.acf[abs(j2)+1]-2*true.acf[j+1]*true.acf[1:lmax+1] W[i,j]<-sum(ti*tj) } } diag(W) @ Therefore, according to this theory, for the lag-$h$ ACF, we have \[ \sqrt{n} \frac{(\widehat \rho_X(h) - \rho_X(h)) }{\sqrt{w_{hh}}} \CiD Normal(0,1) \] which we can assess in our simulation: <>= w.vec.root<-sqrt(diag(W)) std.rho.hat.final<-t(apply(std.rho.hat[,-1],1,function(x){return(x/w.vec.root)})) par(mar=c(4,4,1,2)) boxplot(std.rho.hat.final,col='gray',pch=19,cex=0.7);abline(h=0,lty=2,col='red') apply(std.rho.hat.final,2,var) #Compute the sampling variance of the estimator @ This plot and collection of numerical values demonstrate that the asymptotic theory is holding for this sample size $(n=2000)$ as the sampling variances after standardization are all near one. In a single data set analysis, in a real application, we would have to estimate the $\bW$ matrix also. <>= set.seed(237);n<-2000;lmax<-20 gamma.hat.1<-rho.hat.1<-rep(0,lmax+1) X<-arima.sim(n,model=list(ar=phi.vec,ma=theta.vec)) hmax<-20 W.hat<-matrix(0,hmax,hmax) lmax<-30 est.acf<-c(acf(X,plot=F,lag.max=lmax)$acf,rep(0,1000)) est.acf[1:30] for(i in 1:hmax){ i1<-c(1:lmax)+i i2<-c(1:lmax)-i ti<-est.acf[i1+1]+est.acf[abs(i2)+1]-2*est.acf[i+1]*est.acf[1:lmax+1] for(j in 1:hmax){ j1<-c(1:lmax)+j j2<-c(1:lmax)-j tj<-est.acf[j1+1]+est.acf[abs(j2)+1]-2*est.acf[j+1]*est.acf[1:lmax+1] W.hat[i,j]<-sum(ti*tj) } } diag(W.hat) @ In this calculation, we have assumed that the ACF is zero beyond lag 30; this is because the Bartlett Formula uses ACF values at all positive lags, but the estimation becomes more variable quite quickly. \bigskip The estimation of the autocovariances is, however, much more challenging when the sample size is smaller. Consider the following simulation study when $n=200$; <>= true.acvs<-true.var*ARMAacf(ar=phi.vec,ma=theta.vec,lag.max=20) true.acf<-true.acvs/true.var set.seed(237);n<-200;nits<-5000;lmax<-20 gamma.hat<-rho.hat<-matrix(0,nrow=nits,ncol=lmax+1) for(iter in 1:nits){ X<-arima.sim(n,model=list(ar=phi.vec,ma=theta.vec)) gamma.hat[iter,]<-acf(X,lag.max=lmax,plot=F,type='covariance')$acf rho.hat[iter,]<-gamma.hat[iter,]/gamma.hat[iter,1] } std.gamma.hat<-t(apply(gamma.hat,1,function(x){return(sqrt(n)*(x-true.acvs))})) #Standardize std.rho.hat<-t(apply(rho.hat,1,function(x){return(sqrt(n)*(x-true.acf))})) #Standardize w.vec.root<-sqrt(diag(W)) std.rho.hat.final<-t(apply(std.rho.hat[,-1],1,function(x){return(x/w.vec.root)})) par(mar=c(3,4,0,2));boxplot(std.rho.hat.final,col='gray',pch=19,cex=0.7) abline(h=0,lty=2,col='red') apply(std.rho.hat.final,2,mean) #Expected value of the estimator (should be zero) apply(std.rho.hat.final,2,var) #Sampling variance of the estimator (should be one) @ It is clear from these figures that there is considerable \textbf{bias} in the estimation. \bigskip \textbf{EXAMPLE:} Consider the $\AR(2)$ model \[ X_t = 0.2 X_{t-1} + 0.2 X_{t-2} + Z_t. \] For this model, the polynomial $\Phi(\B)$ is \[ \Phi(\B) = 1 - 0.2 \B - 0.2 \B^2 = (1-\xi_1 \B) (1-\xi_2 \B) \] where $\xi_1 + \xi_2 = 0.2$, $\xi_1 \xi_2 = -0.2$. Note that \[ \Phi(z) = 1 - \phi_1 z - \phi_2 z^2 = 0 \qquad \Longrightarrow \qquad \frac{1}{\xi_1 \xi_2} - \frac{\xi_1 + \xi_2}{\xi_1 \xi_2} z - z^2 = 0 \] provided $|\xi_1|, |\xi_2| > 0$, and hence \[ \left(z-\frac{1}{\xi_1} \right) \left(z - \frac{1}{\xi_2}\right) = 0. \] This we can find the reciprocal roots using \texttt{polyroot} in \texttt{R}, as follows <>= xi<-1/polyroot(c(1,-0.2,-0.2)) @ From this we deduce that $\xi_1$ and $\xi_2$ are $\xi_1 = \Sexpr{round(xi[1],4)} $, $\xi_2 = \Sexpr{round(xi[2],4)}$ and the model is stationary. Consider the following simulation study when $n=200$; <>= phi.vec<-c(0.2,0.2);sigmaZ<-1 true.var<-sigmaZ^2*(1+sum(ARMAtoMA(ar=phi.vec,lag.max=2000)^2)) true.acvs<-true.var*ARMAacf(ar=phi.vec,lag.max=20) true.acf<-true.acvs/true.var set.seed(237);n<-200;nits<-5000;lmax<-20 gamma.hat<-rho.hat<-matrix(0,nrow=nits,ncol=lmax+1) for(iter in 1:nits){ X<-arima.sim(n,model=list(ar=phi.vec)) gamma.hat[iter,]<-acf(X,lag.max=lmax,plot=F,type='covariance')$acf rho.hat[iter,]<-gamma.hat[iter,]/gamma.hat[iter,1] } std.gamma.hat<-t(apply(gamma.hat,1,function(x){return(sqrt(n)*(x-true.acvs))})) #Standardize std.rho.hat<-t(apply(rho.hat,1,function(x){return(sqrt(n)*(x-true.acf))})) #Standardize hmax<-20;W<-matrix(0,hmax,hmax);lmax<-2000 true.acf<-ARMAacf(ar=phi.vec,lag.max=lmax+hmax) for(i in 1:hmax){ i1<-c(1:lmax)+i; i2<-c(1:lmax)-i ti<-true.acf[i1+1]+true.acf[abs(i2)+1]-2*true.acf[i+1]*true.acf[1:lmax+1] for(j in 1:hmax){ j1<-c(1:lmax)+j; j2<-c(1:lmax)-j tj<-true.acf[j1+1]+true.acf[abs(j2)+1]-2*true.acf[j+1]*true.acf[1:lmax+1] W[i,j]<-sum(ti*tj) } } w.vec.root<-sqrt(diag(W)) std.rho.hat.final<-t(apply(std.rho.hat[,-1],1,function(x){return(x/w.vec.root)})) par(mar=c(4,4,1,2)) boxplot(std.rho.hat.final,col='gray',pch=19,cex=0.7);abline(h=0,lty=2,col='red') apply(std.rho.hat.final,2,mean) #Expected value of the estimator (should be zero) apply(std.rho.hat.final,2,var) #Sampling variance of the estimator (should be one) @ Here there is less bias in the estimation, and the asymptotic approximation is somewhat better. \bigskip \textbf{EXAMPLE:} Consider the $\MA(3)$ model \[ X_t = Z_t + 2 Z_{t-1} + 4 Z_{t-2} - 3 Z_{t-3} \] Here the autocorrelation is zero for lags $h > 3$. Consider the following simulation study when $n=200$; <>= theta.vec<-c(2,4,-3);sigmaZ<-1 true.var<-sigmaZ^2*(1+sum(ARMAtoMA(ma=theta.vec,lag.max=2000)^2)) true.acvs<-true.var*ARMAacf(ma=theta.vec,lag.max=20) true.acf<-true.acvs/true.var set.seed(237);n<-200;nits<-5000;lmax<-20 gamma.hat<-rho.hat<-matrix(0,nrow=nits,ncol=lmax+1) for(iter in 1:nits){ X<-arima.sim(n,model=list(ma=theta.vec)) gamma.hat[iter,]<-acf(X,lag.max=lmax,plot=F,type='covariance')$acf rho.hat[iter,]<-gamma.hat[iter,]/gamma.hat[iter,1] } std.gamma.hat<-t(apply(gamma.hat,1,function(x){return(sqrt(n)*(x-true.acvs))})) #Standardize std.rho.hat<-t(apply(rho.hat,1,function(x){return(sqrt(n)*(x-true.acf))})) #Standardize hmax<-20;W<-matrix(0,hmax,hmax) lmax<-2000 true.acf<-ARMAacf(ma=theta.vec,lag.max=lmax+hmax) for(i in 1:hmax){ i1<-c(1:lmax)+i i2<-c(1:lmax)-i ti<-true.acf[i1+1]+true.acf[abs(i2)+1]-2*true.acf[i+1]*true.acf[1:lmax+1] for(j in 1:hmax){ j1<-c(1:lmax)+j j2<-c(1:lmax)-j tj<-true.acf[j1+1]+true.acf[abs(j2)+1]-2*true.acf[j+1]*true.acf[1:lmax+1] W[i,j]<-sum(ti*tj) } } w.vec.root<-sqrt(diag(W)) std.rho.hat.final<-t(apply(std.rho.hat[,-1],1,function(x){return(x/w.vec.root)})) par(mar=c(4,4,1,2));boxplot(std.rho.hat.final,col='gray',pch=19,cex=0.7) abline(h=0,lty=2,col='red') apply(std.rho.hat.final,2,mean) #Expected value of the estimator (should be zero) apply(std.rho.hat.final,2,var) #Sampling variance of the estimator (should be one) @ Here there almost no in the bias in the estimation of the acf, and the asymptotic approximation is somewhat better, although the asymptotic variance seems to be too small as $h$ increases. \end{document}