\documentclass[notitlepage]{article} \usepackage{../Math545} \usepackage{listings} \usepackage{numprint} \usepackage{enumerate} \usepackage{bbm} \usepackage{amsfonts,amsmath} \usepackage{fancyhdr} \usepackage{lastpage} \usepackage{undertilde} \parindent 0in \pagestyle{empty} \setlength{\headheight}{14pt} \def\ExNumber{1} \lfoot{MATH 545 ASSIGNMENT \ExNumber} \cfoot{ } \rfoot{{\small{\em Page \thepage \ of \pageref{LastPage}}}} \renewcommand{\headrulewidth}{0.0pt} \renewcommand{\footrulewidth}{0.0pt} \input{../tcilatex} \def\Var{\text{Var}} \def\Cov{\text{Cov}} \def\Corr{\text{Corr}} \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}} \begin{document} \begin{center} {\textsclarge{MATH 545 - Assignment \ExNumber }} \bigskip {\textbf{ASSIGNMENT}} \end{center} \bigskip In the following questions, $\{Z_t\}\sim \WN(0,\sigma^2)$, with $t \in \mathbb{Z}$, that is, we consider the process defined for all integers \bigskip \begin{enumerate} \item Using the techniques established for the autoregressive model of order 1 (AR(1)), compute the autocorrelation function, $\rho_X(h), h=0,\pm 1, \pm 2, \ldots$ for the following stationary processes: \begin{enumerate} \item $\{X_t\}$ defined implicitly for all $t$ by the equation \[ X_t + 0.1 X_{t-1} = Z_t \] \hfill \textsc{2 Marks} \item $\{X_t\}$ defined implicitly for all $t$ by the equation \[ X_t = 1.1 X_{t-1} - 0.3 X_{t-2} + Z_t \] \textit{Hint: first find the representation \[ X_t = \sum_{j=0}^\infty \psi_j Z_{t-j} \qquad t \in \mathbb{Z} \] for constants $\psi_j$, $j=0,1,2,\ldots$ using the recursive definition. } \hfill \textsc{4 Marks} \item $\{X_t\}$ defined explicitly for all $t$ by the equation \[ X_t = Z_t + Z_{t-1} - 0.5 Z_{t-2} \] \hfill \textsc{3 Marks} \end{enumerate} In each case, show the details of the derivation. \medskip \item The \textit{backshift} operator, $\B$, is an operator that acts on time series processes $\{X_t\}$ such that for every $t$, $\B X_t = X_{t-1}$. The operator has standard algebraic and arithmetic properties such that, for example, we can form polynomials in $\B$: if $\Phi(\B)$ is such a polynomial of order $p \geq 1$, then \[ \Phi(\B) X_t \equiv \left(\sum_{j=0}^p \phi_j \B^j \right) X_t = (\phi_0 + \phi_1 \B + \phi_2 \B^2 + \cdots + \phi_p \B^p) X_t \equiv \sum_{j=0}^p \phi_j X_{t-j}. \] for constants $\phi_0,\phi_1,\ldots,\phi_p$. \medskip Show that the process $\{X_t\}$ from 1 (b) can be written in the form \[ \Phi(\B) X_t = Z_t \] where $\Phi(z) = 1 + \phi_1 z + \phi_2 z^2$ is a polynomial of order two. Show that \[ \Phi(\B) = (1-\alpha_1 \B) (1-\alpha_2 \B) \] for two possibly complex values $\alpha_1,\alpha_2$. \hfill \textsc{5 Marks} \item Using \texttt{R}, generate and plot realizations of length $n=200$ from each of the processes in question 1, assuming that $Z_t \sim Normal(0,1)$. Show the code you use, and provide the resulting plots. \hfill \textsc{6 Marks} \end{enumerate} \pagebreak \bigskip \textit{\textbf{SOLUTIONS}} \begin{enumerate} \item \begin{enumerate} \item We may re-write \[ X_t = - 0.1 X_{t-1} + Z_t \] and recognize this as an $\AR(1)$ process with parameter $\phi = -0.1$. Therefore using the recursion techniques from lectures (using the result that the process is stationary as given in the question) we have that \[ X_t = \sum_{j=0}^\infty \phi^j Z_{t-j} \] and hence the process is zero mean, with \[ \gamma_X(h) = \E[X_t X_{t+h}] = \E \left[\left(\sum_{j=0}^\infty \phi^j Z_{t-j} \right) \left(\sum_{k=0}^\infty \phi^k Z_{t-k} \right) \right] = \E \left[\left(\sum_{j=0}^\infty \phi^j \phi^{j+h} Z_{t-j}^2 \right) \right] = \sigma_Z^2 \phi^h \sum_{j=0}^\infty \phi^{2j} \] as $\E[Z_{t-j} Z_{t+h-k}] = 0$ unless $k = j+h$. Thus, summing the geometric expression, we have \[ \gamma_X(h) = \frac{\sigma_Z^2 \phi^h}{1-\phi^2} \] and hence \[ \rho_X(h) = \frac{\gamma_X(h)}{\gamma_X(0)} = \phi^h = (-0.1)^h. \] This is easily computed in \texttt{R}: the acf for $h=0,1,\ldots,10$ are <>= rho<--0.1 rho^c(0:10) @ \hfill \textsc{2 Marks} \item Several different approaches are possible: \begin{enumerate}[(i)] \item \textbf{Using the infinite sum representation directly:} with the posited form of the solution \[ X_t = \sum_{j=0}^\infty \psi_j Z_{t-j} \] for some constants $\psi_j, j \geq 0$, such that the summation is convergent in mean square (we know that the infinite sum converges in mean square as the process is given as stationary), we can substitute in to the AR formula to obtain \[ \sum_{j=0}^\infty \psi_j Z_{t-j} - \phi_1 \sum_{j=0}^\infty \psi_j Z_{t-j-1} - \phi_2 \sum_{j=0}^\infty \psi_j Z_{t-j-2} = Z_t \] from which we can identify the coefficients, because the result must hold for all $t$. We have \[ \begin{array}{ccc} \text{Term} & \text{LHS} & \text{RHS} \\ Z_t & \psi_0 & 1 \\ Z_{t-1} & \psi_1 - \phi_1 \psi_0 & 0 \\ Z_{t-2} & \psi_2 - \phi_1 \psi_1 - \phi_2 \psi_0 & 0\\ Z_{t-3} & \psi_3 - \phi_1 \psi_2 - \phi_2 \psi_1 & 0\\ \vdots & \vdots\\ Z_{t-j} & \psi_j - \phi_1 \psi_{j-1} - \phi_2 \psi_{j-2} & 0\\ \vdots & \vdots\\ \end{array} \] Thus we have the recursion for $j \geq 2$ \[ \psi_j = \phi_1 \psi_{j-1} + \phi_2 \psi_{j-2} \] with $\psi_0 = 1$ and $\psi_1 = \phi_1$. <>= nlags<-1000; psi.vec<-numeric(nlags) phi1<-1.1 phi2<--0.3 psi0<-1 #psi_0 = 1 psi.vec[1]<-1.1 #psi_1 = 1.1 psi.vec[2]<-phi1*psi.vec[1]+phi2*psi0 for(j in 3:nlags){ psi.vec[j]<-phi1*psi.vec[j-1]+phi2*psi.vec[j-2] #psi_r, r > 1 } psi.vec[1:12] #The first 12 psi values starting at psi_1 ARMAtoMA(ar = c(1.1,-0.3),lag.max=12) #Check using the ARMAtoMA function psi.vec[nlags] #psi_{999} value @ Then using the result (proved from first principles in lectures) that \[ \gamma_X(h) = \sigma_Z^2 \sum_{j=0}^\infty \psi_j \psi_{j+h} \qquad \rho_X(h) = \dfrac{\sum\limits_{j=0}^\infty \psi_j \psi_{j+h}}{\sum\limits_{j=0}^\infty \psi_j^2} \] we can compute the numerical acvs and acfs. Truncating the infinite sum by assuming that $\psi_j = 0$ for $j \geq \Sexpr{nlags}$ yields a computational formula for the covariances, we have: <>= acvs<-numeric(nlags) psi0.vec<-c(psi0,psi.vec) acvs[1]<-sum(psi0.vec^2) #h=0 for(i in 2:nlags){ acvs[i]<-sum(psi0.vec[1:(nlags-i+1)]*psi0.vec[i:nlags]) #h=i } acfs<-acvs/acvs[1] acvs[1:10] #The first 10 acvs values acfs[1:10] #The first 10 acfs values @ Note that the recursion above can be solved using linear recurrence relation methods: we have \[ \left[ \begin{array}{c} \psi_j \\ \psi_{j-1} \end{array} \right] = \left[ \begin{array}{cc} \phi_1 & \phi_2 \\[4pt] 1 & 0 \\[4pt] \end{array} \right] \left[ \begin{array}{c} \psi_{j-1} \\ \psi_{j-2} \end{array} \right] = \left[ \begin{array}{cc} \phi_1 & \phi_2 \\[4pt] 1 & 0 \\[4pt] \end{array} \right]^{j-1} \left[ \begin{array}{c} \psi_{1} \\ \psi_{0} \end{array} \right] \] <>= nlags<-1000; psi.vec<-numeric(nlags) Phi<-matrix(c(phi1,1,phi2,0),2,2) psi.mat<-matrix(0,ncol=nlags,nrow=2) psi.mat[,1]<-c(1.1,1) #initialize the recursion for(j in 2:nlags){ psi.mat[,j]<-Phi %*% psi.mat[,j-1] #apply the recursive formula } psi.mat[1,1:12] @ \item \textbf{Using the recurrence relation:} we want to obtain $X_t$ explicitly in terms of $\{Z_t\}$. We may use the recursion method to compute such a representation: to begin, we substitute in the equivalent expression for $X_{t-1}$ \begin{align} X_t = \phi_1 X_{t-1} + \phi_2 X_{t-2} + Z_t & = \phi_1 (\phi_1 X_{t-2} + \phi_2 X_{t-3} + Z_{t-1}) + \phi_2 X_{t-2} + Z_t \nonumber \\[6pt] & = \alpha_1^{(2)} X_{t-2} + \alpha_2^{(2)} X_{t-3} + \sum_{j=0}^1 \psi_j Z_{t-j} \label{eq:step1} \end{align} say, where $\psi_1 = 0, \psi_2 = \phi_1$ and \[ \alpha_1^{(2)} = \phi_1^2 + \phi_2 \qquad \alpha_2^{(2)} = \phi_1 \phi_2 \] Let $\mathbf{\alpha}^{(2)} = (\alpha_1^{(2)},\alpha_2^{(2)},\ldots,\alpha_p^{(2)})^\top$ be the vector of coefficients. Then \[ \alpha_1^{(2)} = \phi_1 \alpha_1^{(1)} + \alpha_2^{(1)} \qquad \alpha_2^{(2)} = \phi_2 \alpha_1^{(1)}. \] Therefore we may write $\mathbf{\alpha}^{(2)} = \bA \mathbf{\alpha}^{(1)}$ for $(2 \times 2)$ matrix $\bA$, where \[ \bA = \left[ \begin{array}{cc} \phi_1 & 1 \\[4pt] \phi_2 & 0 \\[4pt] \end{array} \right] \] By repeatedly applying the recursion, we obtain \[ \mathbf{\alpha}^{(r)} = \bA \mathbf{\alpha}^{(r-1)} = \bA^2 \mathbf{\alpha}^{(r-2)} = \cdots = \bA^{r-1} \mathbf{\alpha}^{(1)}. \] which defines the coefficients in the representation \begin{align} X_t & = \alpha_1^{(r)} X_{t-r} + \alpha_2^{(r)} X_{t-r-1} + \sum_{j=0}^{r-1} \psi_j Z_{t-j} \label{eq:step1} \end{align} \[ \mathbf{\alpha}^{(r)} = \bA \mathbf{\alpha}^{(r-1)} \] with $\psi_j = \alpha_1^{(j)}$. Diagonalizing $\bA$, that is, writing \[ \bA = \bV \bD \bV^{-1} \qquad \Longrightarrow \qquad \bA^r = \bV \bD^r \bV^{-1} \] Now the eigenvalues can be found using the characteristic equation \[ \left|\left[ \begin{array}{cc} \phi_1-\lambda & 1 \\[4pt] \phi_2 & -\lambda \\[4pt] \end{array} \right] \right| = \lambda(\lambda-\phi_1) - \phi_2 = 0 \] yielding the solutions $\left(\phi_1 \pm \sqrt{\phi_1^2 + 4 \phi_2}\right)/2$. Numerically, this yields the two eigenvalues $\lambda_1 = 0.6, \lambda_2=0.5$, which are both less than one in modulus, so we have that $\bA \longrightarrow \mathbf{0}_{p,p}$ as required. Once we have computed \[ \mathbf{\alpha}^{(r)} = \bA^{r-1} \mathbf{\alpha}^{(1)} \] for each $r$, we can generate the $\psi_j$ values by recalling that $\psi_j = \alpha_1^{(j)}$, and then using the formula above for the acvs and acfs. Using \texttt{R} for the numerical solution, with matrix $\bA$ \[ \bA = \left[ \begin{array}{rc} 1.1 & 1 \\[4pt] -0.3 & 0 \\[4pt] \end{array} \right] \] we can check the eigen decomposition in \texttt{R} <>= A<-matrix(c(1.1,-0.3,1,0),2,2) eigen(A) @ indicating that the eigenvalues are 0.6, 0.5. Thus, $\bA^r$ should converge to the zero matrix as $r \longrightarrow \infty$; we can check that also <>= A.r<-diag(c(1,1)) #Identity matrix for(r in 1:100){A.r<-A.r %*% A} #A.r stores the A^r matrix A.r @ This demonstrates that the process can be represented \[ X_t = \sum_{j=0}^\infty \psi_j Z_{t-j} \] and that this process is well-defined (in mean-square terms). To compute the $\psi$ values, we may use the recursive formula. We perform the calculation out to 1000 steps: <>= nlags<-1000; psi.vec<-numeric(nlags) psi0<-1 #psi_0 = 1 psi.vec[1]<-1.1 #psi_1 = 1.1 alpha.mat<-matrix(0,ncol=nlags,nrow=2) alpha.mat[,1]<-c(1.1,-0.3) #initialize the recursion for(r in 2:nlags){ alpha.mat[,r]<-A %*% alpha.mat[,r-1] #apply the recursive formula psi.vec[r]<-alpha.mat[1,r] #psi_r, r > 1 } psi.vec[1:12] #The first 12 psi values starting at psi_1 ARMAtoMA(ar = c(1.1,-0.3),lag.max=12) #Check using the ARMAtoMA function psi.vec[nlags] #psi_{999} value @ Once we have the representation in terms of $\psi$ values, we may compute from first principles the autocovariances. <>= acvs<-numeric(nlags) psi0.vec<-c(psi0,psi.vec) acvs[1]<-sum(psi0.vec^2) #h=0 for(i in 2:nlags){ acvs[i]<-sum(psi0.vec[1:(nlags-i+1)]*psi0.vec[i:nlags]) #h=i } acfs<-acvs/acvs[1] acvs[1:10] #The first 10 acvs values acfs[1:10] #The first 10 acfs values @ We can check this result using the \texttt{ARMAacf} function in \texttt{R} <>= ARMAacf(ar=c(1.1,-0.3),lag.max=9) @ Note that an analytic (closed-form) solution to the recurrence relation is possible using the eigenvector matrix $\bV$ as we may write \[ \psi_j = \alpha_1^{(j)} = c_1 \lambda_1^j + c_2 \lambda_2^j \qquad j \geq 0 \] where $c_1,c_2$ are derived from the diagonalization result $\bA^r = \bV \bD^r \bV^{-1}$. \item \textbf{Recursion on the acv:} using the technique from lectures (multiplying both sides by $X_{t-h}$ and taking expectations), we have \[ \E[X_{t} X_{t-h}] = \phi_1 \E[X_{t-1} X_{t-h}] + \phi_2 \E[X_{t-2} X_{t-h}] + \E[Z_t X_{t-h}] \] with $\phi_1 = 1.1$ and $\phi_2 = -0.3$. From this we conclude that \[ \gamma_X(h) = \phi_1 \gamma_X(h-1) + \phi_2 \gamma_X(h-2) \qquad h \geq 2. \] as $Z_t$ and $X_{t-h}$ are uncorrelated by standard results. When $h=0$, we have \[ \gamma_X(0) = \phi_1 \gamma_X(-1) + \phi_2 \gamma_X(-2) \equiv \phi_1 \gamma_X(1) + \phi_2 \gamma_X(2) \] as we know $\gamma_X(-h) = \gamma_X(h)$, and when $h=1$ \[ \gamma_X(1) = \phi_1 \gamma_X(0) + \phi_2 \gamma_X(-1) \equiv \phi_1 \gamma_X(0) + \phi_2 \gamma_X(1). \] From the latter equation we have \[ \gamma_X(1) = \frac{\phi_1}{1-\phi_2}\gamma_X(0). \] Therefore, for each $h \geq 2$, for the acfs \[ \rho_X(h) = \phi_1 \rho_X(h-1) + \phi_2 \rho_X(h-2) \] with $\rho_X(0) = 1, \rho_X(1) = \phi_1/(1-\phi_2)$. This recursion allows ready computation of the acfs. <>= acfs<-numeric(nlags) acfs[1]<-1 #h=0 acfs[2]<-phi1/(1-phi2) #h=1 for(i in 3:nlags){ acfs[i]<- phi1*acfs[i-1] + phi2*acfs[i-2] #h=i } acfs[1:10] #The first 10 acfs values @ Note that, once again, for an analytic solution, we have the linear recurrence relation \[ \left[ \begin{array}{c} \rho_X(h) \\ \rho_X(h-1) \end{array} \right] = \left[ \begin{array}{cc} \phi_1 & \phi_2 \\[4pt] 1 & 0 \\[4pt] \end{array} \right] \left[ \begin{array}{c} \rho_X({h-1}) \\ \rho_X({h-2}) \end{array} \right] = \left[ \begin{array}{cc} \phi_1 & \phi_2 \\[4pt] 1 & 0 \\[4pt] \end{array} \right]^{h-1} \left[ \begin{array}{c} \rho_X({1}) \\ \rho_X({0}) \end{array} \right] \] and the same form of solution as above. \end{enumerate} \hfill \textsc{4 Marks} \item In this case we have the $\MA(2)$ process with coefficients $\theta_0 = \theta_1 = 1$, $\theta_2 = -0.5$. By direct calculation (or citing the result from lectures), we know that $\rho_X(h) = 0$ if $|h| > q = 2$. For the other three cases, we have \[ \gamma_X(h) = \left\{ \begin{array}{cc} \theta_0^2 + \theta_1^2 + \theta_2^2 & h = 0 \\ \theta_0 \theta_1 + \theta_1 \theta_2 & h = \pm 1 \\ \theta_0 \theta_2 & h = \pm 2 \\ \end{array} \right. \] <>= th0<-th1<-1 th2<--0.5 acvs<-numeric(4) acvs[1]<-th0^2 + th1^2 + th2^2 #h=0 acvs[2]<-th0*th1 + th1*th2 #h=1 acvs[3]<-th0*th2 #h=2 acvs acfs<-acvs/acvs[1] acfs ARMAacf(ma=c(1,-0.5),lag.max=3) @ Numerically, the three values are $\rho_X(0) = \Sexpr{acfs[1]}$, $\rho_X(1) = \Sexpr{acfs[2]}$ and $\rho_X(2) = \Sexpr{acfs[3]}$. \hfill \textsc{3 Marks} \end{enumerate} In each case, show the details of the derivation. \medskip \item By rearrangement of the implicit equation, we have that \[ X_t - 1.1 X_{t-1} + 0.3 X_{t-2} = Z_t \] hence $\Phi(\B) = 1 - 1.1 \B + 0.3 \B^2$. Then we see that in fact \[ \Phi(\B) = (1-0.6 \B) (1- 0.5 \B) \] so $\alpha_1=0.6,\alpha_2=0.5$. Any quadratic polynomial can be factored into the product of two factors in this way, but here the $\alpha$ quantities are real-valued. Note that these two values match the eigenvalues from question 1b. \hfill \textsc{5 Marks} \item For the three models, any sensible starting values can be chosen in the $\AR(p)$ case for 1 (a) and (b); for example, it is permissible here to set the initial values to zero (the expected value of $X_t$) or to randomly chosen $Z\sim Normal(0,1)$ values. It is also possible (and in fact preferable) to generate $X_t$ for $t \geq -20$ say, after initializing as mentioned. For (c), we can generate $Z_{-2},Z_{-1},Z_{0}$ and generate the $X_t$ values accordingly. <>= set.seed(2101) n<-200;Z<-rnorm(n+20);X<-rep(0,n+20) phi1<--0.1 X[1]<-0 for(i in 2:(n+20)){ X[i]<-phi1*X[i-1]+Z[i] } X<-X[21:(n+20)] par(mar=c(4,4,0,0)) plot(1:n,X,xlab='t',type='l',ylim=range(-5,5)) points(1:n,X,pch=19,cex=0.5) @ <>= set.seed(2101) n<-220;Z<-rnorm(n+20);X<-rep(0,n+20) phi1<-1.1;phi2<--0.3 for(i in 3:(n+20)){ X[i]<-phi1*X[i-1]+phi2*X[i-2]+Z[i] } X<-X[21:(n+20)] par(mar=c(4,4,0,0)) plot(1:n,X,xlab='t',type='l',ylim=range(-5,5)) points(1:n,X,pch=19,cex=0.5) @ <>= set.seed(2101) #Set the random number generator seed value n<-200 #Set the sample size Z<-rnorm(n) #Generate Z[1],...Z[n] Z.0<-rnorm(1) #Generate Z[0] Z.1<-rnorm(1) #Generate Z[-1] Z.extra<-c(Z.1,Z.0,Z) #Concatenate X<-Z.extra[3:(n+2)]+Z.extra[2:(n+1)]-0.5*Z.extra[1:n] #Define the time series X par(mar=c(4,4,0,0)) plot(1:n,X,xlab='t',type='l',ylim=range(-5,5)) points(1:n,X,pch=19,cex=0.5) @ \hfill \textsc{6 Marks} Note: In this case, as the white noise process is Gaussian, we can in fact proceed using the fact that $\{X_t\}$ must also be Gaussian, with mean zero and variance covariance matrix given by the acvs. Therefore we can use the \texttt{mvrnorm} function in \texttt{R}: we also use the computational functions \texttt{ARMAtoMA} to obtain the linear process expansion and compute the variance of $\{X_t\}$, \texttt{ARMAacf} to compute the acf, and \texttt{toeplitz} to form the matrix. <>= set.seed(2101) library(MASS) n<-200; sigX1<-sum(ARMAtoMA(c(-0.1),lag.max=100)^2) Sig1<-toeplitz(sigX1*ARMAacf(ar=c(-0.1), lag.max = 199)) X<-mvrnorm(n,mu=rep(0,n),Sigma=Sig1) sigX2<-sum(ARMAtoMA(c(1.1,-0.3),lag.max=100)^2) Sig2<-toeplitz(sigX2*ARMAacf(ar=c(1.1,-0.3), lag.max = 199)) X<-mvrnorm(n,mu=rep(0,n),Sigma=Sig2) sigX3<-sum(c(1,1,-0.5)^2) Sig3<-toeplitz(sigX3*ARMAacf(ma=c(1,-0.5), lag.max = 199)) X<-mvrnorm(n,mu=rep(0,n),Sigma=Sig3) @ \end{enumerate} \end{document}