\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}} \newcommand{\bs}[1]{\boldsymbol{\bm{\unslant{#1}}}} \newcommand{\Rho}{\mathrm{P}} \newcommand{\cif}{\text{ if }} \newcommand{\Ra}{\Longrightarrow} \newcommand{\ra}{\longrightarrow} %\def\phivec{\boldsymbol{\bm{\unslant \phi}}} %\def\thetavec{\boldsymbol{\bm{\unslant \theta}}} %\def\rhovec{\boldsymbol{\bm{\unslant \rho}}} %\def\gammavec{\boldsymbol{\bm{\gamma}}} \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{Algorithms for Forecasting} \begin{document} \maketitle Recall that the optimal linear forecast for a stationary process is found by solving a least-squares problem: we minimize the mean square (forecast or prediction) error \begin{equation}\label{eq:MSE1} \MSE_n(h) = \E\left[\left(X_{n+h} - \widehat X_{n+h} \right)^2 \right] \end{equation} for a specified $n$ and $h \geq 1$. We first note the decomposition \begin{align} \E\left[\left(X_{n+h} - \widehat X_{n+h} \right)^2 \right] & = \E\left[\left(X_{n+h} - m_n(h) + m_n(h) - \widehat X_{n+h} \right)^2 \right] \nonumber\\[6pt] & = \E\left[(X_{n+h} - m_n(h))^2 \right] + 2 \E\left[ (X_{n+h} - m_n(h))(m_n(h) - \widehat X_{n+h}) \right] + \E\left[ (m_n(h) - \widehat X_{n+h})^2 \right] \label{eq:MSE2} \end{align} where $m_n(h) = \E[X_{n+h}|X_1,\ldots,X_n]$ is the conditional expectation of $X_{n+h}$ given $X_1,\ldots,X_n$. In these three terms, the first is some constant that does not depend on the choice of $\widehat X_{n+h}$; the second is zero, as, by iterated expectation \[ \E\left[ (X_{n+h} - m_n(h))(m_n(h) - \widehat X_{n+h}) \right] = \E_{X_{1:n}}\left[ \E_{X_{n+h}|X_{1:n}} \left[ (X_{n+h} - m_n(h))(m_n(h) - \widehat X_{n+h}) |X_{1:n} \right] \right] \] if $X_{1:n} \equiv (X_1,\ldots,X_n)$ -- for the inner expectation, we have \[ \E_{X_{n+h}|X_{1:n}} \left[ (X_{n+h} - m_n(h))(m_n(h) - \widehat X_{n+h}) |X_{1:n} \right] = (m_n(h) - \widehat X_{n+h}) \E_{X_{n+h}|X_{1:n}} \left[ (X_{n+h} - m_n(h)) |X_{1:n} \right] \] as both $m_n(h)$ and $ \widehat X_{n+h}$ are functions of $X_{1:n}$, and $\E_{X_{n+h}|X_{1:n}} \left[ (X_{n+h} - m_n(h)) |X_{1:n} \right] = 0$. Therefore we mininmize \eqref{eq:MSE1} by minimizing the third term in \eqref{eq:MSE2}: this can be achieved by setting \[ \widehat X_{n+h} = m_n(h) = \E[X_{n+h}|X_1,\ldots,X_n] \] in which case \begin{equation}\label{eq:MinMSE} \MSE_n(h) = \E\left[\left(X_{n+h} - \E[X_{n+h}|X_1,\ldots,X_n] \right)^2 \right]. \end{equation} Using a linear representation for $m_n(h)$, we can get a more explicit result \[ \widehat X_{n+h} = a_0 + \sum_{i=1}^n a_i X_{n-i+1}. \] This reduces to choosing the coefficients $a = (a_1,\ldots,a_n)^\top$. We minimize using calculus: we have \begin{align} \dfrac{\partial \MSE_n(h)}{\partial a_0} & = \E\left[\left(X_{n+h} - a_0 - \sum_{i=1}^n a_i X_{n-i+1} \right) \right] \label{eq:OLS1}\\[6pt] \dfrac{\partial \MSE_n(h)}{\partial a_j} & = \E\left[\left(X_{n+h} - a_0 - \sum_{i=1}^n a_i X_{n-i+1} \right) X_{n-j+1} \right], j=1,\ldots,n. \label{eq:OLS2} \end{align} allowing for differentiation under the expectation. We must solve after setting \eqref{eq:OLS1} and \eqref{eq:OLS2} equal to zero simultaneously. From equation \eqref{eq:OLS1}, we have that \[ \E[X_{n+h}] - a_0 - \sum_{i=1}^n a_i \E[X_{n-i+1}] = 0 \] and if $\{X_t\}$ is stationary with mean $\mu$, we may conclude that the optimal choice must satisfy \[ a_0 = \mu \left(1-\sum_{i=1}^n a_i \right). \] From \eqref{eq:OLS2}, we have for each $j=1,\ldots,n$, \[ \E[X_{n+h} X_{n-j+1}] - a_0 \E[X_{n-j+1}] - \sum_{i=1}^n a_i \E[X_{n-i+1} X_{n-j+1}] = 0. \] Noting that from the earlier result, we have \[ a_0 \E[X_{n-j+1}] = \mu^2 \left(1-\sum_{i=1}^n a_i \right) \] and hence, rearranging \[ \E[X_{n+h} X_{n-j+1}] - \mu^2 - \sum_{i=1}^n a_i (\E[X_{n-i+1} X_{n-j+1}] - \mu^2) = 0. \] so that \[ \sum_{i=1}^n a_i \Cov[X_{n-i+1},X_{n-j+1}] = \Cov[X_{n+h},X_{n-j+1}] \] or equivalently \[ \sum_{i=1}^n a_i \gamma_X(i-j) = \gamma_X(j+h-1)\qquad j=1,\ldots,n. \] Stacking these $n$ equations, we are left with \[ \begin{bmatrix} \gamma_X(0) & \gamma_X(1) & \gamma_X(2) & \cdots & \gamma_X(n-1) \\[6pt] \gamma_X(1) & \gamma_X(0) & \gamma_X(1) & \cdots & \gamma_X(n-2)\\[6pt] \gamma_X(2) & \gamma_X(1) & \gamma_X(0) & \cdots & \gamma_X(n-2)\\[6pt] \vdots & \vdots & \vdots & \ddots & \vdots \\[6pt] \gamma_X(n-1) & \gamma_X(n-2) & \gamma_X(n-3) & \cdots & \gamma_X(0) \end{bmatrix} \begin{bmatrix} a_1 \\[6pt] a_2 \\[6pt] a_3 \\ \vdots \\[6pt] a_n \end{bmatrix} = \begin{bmatrix} \gamma_X(h)\\[6pt] \gamma_X(h+1) \\[6pt] \gamma_X(h+2)\\[6pt] \vdots \\[6pt] \gamma_X(h+n-1) \end{bmatrix} \] or equivalently \begin{equation}\label{eq:LinSys} \Gammat_X(n) \bs{a} = \Gammat_X(n,h) \end{equation} say, where \[ \Gammat_X(n,h) = (\gamma_X(h),\ldots,\gamma_X(h+n-1))^\top. \] Thus provided $\Gammat_X(n)$ is non-singular, we have that the optimal $\bs{a}$ is the solution to \eqref{eq:LinSys} \[ \bs{a}(n,h) = \{ \Gammat_X(n) \}^{-1} \gammavec_X(n,h) \] which requires us to compute $\{ \Gammat_X(n) \}^{-1}$. \pagebreak \textbf{The Levinson-Durbin Algorithm:} Without loss of generality, take $\E[X_t]=0$. Consider forecasting for $h=1$, and write \[ \widehat{X}_{n+1}=\sum_{i=1}^n \phi_{n,i}X_{n-i+1} = \phivec_n^\top X_{n:1} \] where $\phivec_n=(\phi_{n,1}, \cdots, \phi_{n,n})^\top$, which corresponds in our prediction forecast to $\bs{a}(n,1)$. The optimal choice solves the linear system \begin{align*} \Gammat_X(n)\phivec_n=\Gammat_X(n,1) = \gammavec_n^\tR, \end{align*} say, which is an $n \times 1 $ system of equations, yielding minimum MSE \begin{align*} v_n=\gamma_X(0)-\phivec_n^\top \gammavec_n^\tR \end{align*} The Levinson-Durbin recursion solves the equation without matrix inversion, in order $n^2$ operations (most matrix inversion procedures are order $n^3$) by exploiting the Toeplitz structure of $\Gamma_X(n)$. One can also use singular value, eigenvalue decomposition or Cholesky decomposition; these method have special implementation for positive-definite Toeplitz structure matrices. \bigskip \emph{Initialization}: For $n=1$, set \[ \phi_{1,1}=\frac{\gamma_X(1)}{\gamma_X(0)}=\rho_X(1) \qquad v_0=\gamma_X(0) \qquad v_1=\gamma_X(0)(1-\phi_{1,1}^2). \] \medskip \emph{Recursion:} Define $\phivec_n$ in terms of $\phivec_{n-1}$ and $v_{n-1}$ as follows: \begin{align*} \phi_{n,n}= \frac{1}{v_{n-1}} \left[ \gamma_X(n) - \phivec^\top_{n-1} \gammavec_{n-1}^\tRb \right] \end{align*} where $\Gammat_X^\tRb(n-1,1)=(\gamma_X(n-1), \gamma_X(n-1), \ldots, \gamma_X(1))^\top $ (the superscript `R' means `reverse') and \[\phivec_{n,1:(n-1)}=\phivec_{n-1}-\phi_{n,n} \phivec_{n-1}^{\tRb}\] where $\phivec_{n-1}$ and $\phivec_{n-1}^{\tRb}$ are $(n-1)\times 1$ vectors with $\phivec_{n-1}^{\tRb}=(\phi_{n-1,n-1}, \phi_{n-1, n-2}, \cdots,\phi_{n-1,1})^\top$, and $\phi_{n,n}$ is a scalar Then set \[ v_n=v_{n-1}(1-\phi^2_{n,n}). \] \medskip \begin{proof} \smallskip To see why this works, let \begin{align*} \bR_n=\frac{1}{\gamma_X(0)}\Gammat_X(n) \end{align*} be the autocorrelation matrix and let \[ \bs{\rho}_{n}= \frac{1}{\gamma_X(0)} (\gamma_X(1), \ldots, \gamma_X(n))^\top \qquad \bs{\rho}_{n}^\tRb = (\rho_X(n), \cdots, \rho_X(1))^\top. \] The recursion computes the inverse of $\bP_{n+1}$ from the inverse of $\bP_n$ using rank-one operations. Note first that \[ \Gammat_X(n) \phivec_n=\gammavec_n^\tRb \Ra\bR_n \phivec_n = \bs{\rho}_{n}. \] We need to verify that \[ \bR_n \phivec_n = \bs{\rho}_{n} \Ra \bR_{n+1}\phivec_{n+1}=\bs{\rho}_{n+1}. \] For $n=1$, $\bR_1=1,\phivec_{1}=\phi_{1,1}=\rho(1)=\bs{\rho}_{1,1}$. We establish the result by induction on $n.$ Assume that $\bR_k\phivec_k=\bs{\rho}_{k}$. Now \begin{align*} \bR_{k+1}=\begin{bmatrix} \bR_k &\bs{\rho}_{k}^\tR\\[6pt] \{\bs{\rho}_{k}^\tR\}^\top &1 \end{bmatrix} \end{align*} and $\bR_k\bs{\phi_k}=\bs{\rho}_{k}$, therefore \begin{align*} \bR_k \phivec_k^{\tRb}=\bs{\rho}_{k}^\tRb \end{align*} using the symmetry of $\bR_k$. In the proposed recursion, we have that \begin{align*} \bR_{k+1}\phivec_{k+1}= \begin{bmatrix} \bR_k &\bs{\rho}_{k}^\tR\\[6pt] \{\bs{\rho}_{k}^\tR\}^\top &1 \end{bmatrix} \begin{bmatrix} \phivec_k-\phi_{k+1,k+1}\phivec_{k}^{\tRb} \\[6pt] \phi_{k+1,k+1} \end{bmatrix} = \begin{bmatrix} \ba\\[6pt] \bb \end{bmatrix} \end{align*} where $ \ba$ is equal to a $k\times 1$ system of the form \[ \ba=\bR_k\phivec_k - \phi_{k+1,k+1}\bR_k\phivec_k^{\tRb} + \phi_{k+1,k+1}\bs{\rho}_{k}=\bR_k\phivec_k=\bs{\rho}_{k} \] by the inductive hypothesis. For $ \bb$, we have \[ \bb=\{\bs{\rho}_{k}^\tR\}^\top \phivec_k - \phi_{k+1,k+1} \{\bs{\rho}_{k}^\tR\}^\top \phivec_k^{\tRb} + \phi_{k+1, k+1} = \{\bs{\rho}_{k}^\tR\}^\top \phivec_k - \phi_{k+1,k+1}(1-\{\bs{\rho}_{k}^\tR\}^\top \phivec_k^{\tRb}) \] %Equal to $\rho(k+1)$ We have now \[ \phi_{k+1,k+1}= \frac{1}{v_k} \left[\gamma_X(k+1)-\phivec_k^\top \gammavec_k^{\tRb} \right] \] where \[ v_k =\gamma_X(0)-\phivec_k^\top \Gammat_X(k,1)=\gamma_X(0) (1-\phivec_k^\top \bs{\rho}_k ) = \gamma_X(0) (1- \{\bs{\rho}_{k}^\tRb\}^\top \phivec_k^{\tRb}) \] From the above formula, we have \begin{align} \phi_{k+1,k+1}=\frac{1}{\gamma_X(0) (1-\{\bs{\rho}_{k}^\tR\}^\top \phivec_k^{\tR})} \left[\gamma_X(k+1)-\phivec_k^\top \gammavec_k^\tRb \right] \label{Levinson-D} \end{align} changing the order of the inner product and from $\mathbf{b}$, we have \begin{align*} \bs{\rho}_k^\top \phivec_k + \frac{\left[\gamma_X(k+1)-\phivec_k^\top \gammavec_k^\tRb \right]}{\gamma_X(0) (1-\bs{\rho}_k^\top \phivec_k^{\tR})} (1-\bs{\rho}^\top_{k} \phivec_k^{ \tRb }) &= \{\bs{\rho}_{k}^\tRb\}^\top \phivec_k +\rho_X(k+1)-\phivec_k^\top \bs{\rho}_k^\tRb\\[6pt] &=\rho_X(k+1) \end{align*} as required. Therefore the required recursion is valid by the inductive argument. \end{proof} \medskip We have effectively performed a matrix inversion by considering a block decomposition \begin{align*} \bR_{k+1}=\begin{bmatrix} \bR_k &\bs{\rho}_{k}^\tR\\[6pt] \{\bs{\rho}_{k}^\tR\}^\top &1 \end{bmatrix} \end{align*} and then computed $\bR_{k+1}^{-1}$ using $\bR_k^{-1}$. This works because of the nature of the nature of the matrix (Toeplitz structure). \medskip For $v_n$, we have \begin{align*} v_n=\E[(X_{n+1}-\phivec_n^\top X_{n:1})^2]&=\gamma_X(0)-\phivec_n^\top \gammavec_{n-1}^\tRb \\ & = \gamma_X(0)-\phivec_{n-1}^\top \gammavec_{n-1}^\tRb +\phi_{n,n} \{\phivec_{n-1}^\tRb \}^\top \gammavec_{n-1}^\tRb-\phi_{n,n}\gamma_X(n)\\ & = v_{n-1}+\phi_{n,n}\big( \{\phivec_{n-1}^\tRb\}^\top \gammavec_{n-1}^\tRb -\gamma_X(n)\big) \\&=v_{n-1}-\phi_{n,n}^2 \big(\gamma_X(0)- \{\phivec_{n-1}^\tRb\}^\top \gammavec_{n-1}^\tRb \big) \\&= v_{n-1}-\phi_{n,n}^2 v_{n-1} \\&= v_{n-1}(1-\phi^2_{n,n}) \end{align*} plugging in \eqref{Levinson-D}. That is, $v_n=v_{n-1}(1-\phi_{n,n}^2)$ and as $v_n \geq 0$, imply $\phi_{n,n}^2<1, |\phi_{n,n}|<1$ thus $v_n \leq v_{n-1}$, that is, the $\{v_n\}$ form a non-increasing sequence. \pagebreak \textbf{Partial autocorrelation:} Using the Levinson-Durbin algorithm, define the function $\alpha_X(h)$ by \begin{align*} \alpha_X(h)= \begin{cases} 1& \cif h=0\\\phi_{h,h} & \cif h \geq 1 \end{cases} \end{align*} for $h=0,1,2,\ldots $ This $\alpha_X(h)$ is the \emph{partial autocorrelation function} (PACF). For random variables $X,Y,Z$, the \textit{partial correlation} is defined by \begin{align*} \Corr [(X-\E[X|Z]),(Y-\E[Y|Z])]; \end{align*} this is the correlation between residuals obtained by regressing in turn $X$ and on $Z$ and $Y$ on $Z$. The PACF computes the autocorrelation between \[ X_t-\E[X_t|X_{(t+1):(t+h-1)}] \qquad \text{and} \qquad X_{t+h}-\E[X_{t+h}|X_{(t+1):(t+h-1)}]. \] \medskip Consider a causal stationary AR$(p)$ process; for $h>p$, \[ \E[X_{t+h}|X_{(t+1):(t+h-1)}]=\phi_1 X_{t+h-1}+\cdots + \phi_pX_{t+h-p} \qquad \E[X_{t}|X_{(t+1):(t+h-1)}]= f(X_{t+1}, \ldots , X_{t+h-1}) \] for some function $f()$, where by definition $X_{t+h}-\E[X_{t+h}|X_{(t+1):(t+h-1)}]=Z_{t+h}$, therefore \[ \Corr [X_t-\E[X_t|X_{(t+1):(t+h-1)}],X_{t+h}-\E[X_{t+h}|X_{(t+1):(t+h-1)}]] = \Corr [X_t-\E[X_t|X_{(t+1):(t+h-1)}],Z_{t+h}] = 0 \] as $\{X_t\}$ is causal so $X_t-\E[X_t|X_{(t+1):(t+h-1)}]$ depends on $\{Z_s, s \leq t+h-1\}$ only. Thus $\alpha_X(h)=0$ for $h>p$, while for $h \leq p, \alpha_X(h) \neq 0. $ That is, we can diagnose an $\AR(p)$ structure by inspecting whether the partial autocorrelation drops to zero at some finite lag. \bigskip \textbf{EXAMPLE} Consider the $\AR(3)$ process \[ X_t = -1.4 X_{t-1} - 0.48 X_{t-2} + 0.16 X_{t-3} + Z_t \] or equivalently \[ (1 + 1.4 \B + 0.48 \B^2 - 0.16 \B^3) X_t = Z_t. \] or $\Phi(B)X_t = Z_t$ with $\Phi(z) = (1 + 1.4 z + 0.48 z^2 - 0.16 z^3)$. We have the factorization \[ \Phi(z) = (1-0.2z)(1+1.6z+0.8z^2) = (1-\xi_1z)(1-\xi_2 z) (1-\xi_3 z) \] and first reciprocal root is $\xi_1 = 0.2$. The two remaining roots are \[ \frac{-1.6 \pm \sqrt{2.56 - 3.2}}{2} = \frac{-1.6 \pm \sqrt{-0.64}}{2} = -0.8 \pm 0.4 \imagi \] so the reciprocal roots are $\xi_2 = -0.8 + 0.4 \imagi$ and $\xi_3 = -0.8 - 0.4\imagi$. All the reciprocal roots are less than one in modulus, so the process is stationary and causal. <>= alpha.vals<-c(0.2,-0.8+0.4i,-0.8-0.4i) 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] } ar.coefs #These are the AR coefficients n<-200 rho.fn1<-ARMAacf(ar=Re(-ar.coefs[2:4]), lag.max = n+1) #Simulate data par(mar=c(4,4,0,1)) plot(rho.fn1[1:50],type="n",ylim=range(-1,1),xlab="h",ylab=expression(rho[X](h))) abline(h=0,lty=2) for(h in 1:50){lines(c(h,h),c(0,rho.fn1[h]))} points(rho.fn1[1:50],pch=19,cex=0.5) @ <>= set.seed(26) y<-arima.sim(n=1000,model=list(ar=Re(-ar.coefs[2:4]))) par(mar=c(4,4,0,1)) plot(y) ################################################## #Levinson-Durbin psi.Vals<-ARMAtoMA(ar=Re(-ar.coefs[2:4]),lag.max=10000) Variance<-1+sum(psi.Vals^2) gamma.fn1<-Variance*rho.fn1[-1] Phi.mat<-matrix(0,nrow=n,ncol=n) v<-rep(0,n+1) v[0+1]<-Variance Phi.mat[1,1]<-gamma.fn1[1]/Variance #Initialize v[1+1]<-Variance*(1-Phi.mat[1,1]^2) for(j in 2:n){ #Recursion Phi.mat[j,j]<-(gamma.fn1[j]-sum(Phi.mat[j-1,1:(j-1)]*rev(gamma.fn1[1:(j-1)])))/v[(j-1)+1] Phi.mat[j,1:(j-1)]<- Phi.mat[j-1,1:(j-1)] - Phi.mat[j,j] * rev(Phi.mat[j-1,1:(j-1)]) v[j+1]<-v[(j-1)+1]*(1-Phi.mat[j,j]^2) } round(Phi.mat[1:8,1:8],10) diag(Phi.mat)[1:10] #pacfs ARMAacf(ar=Re(-ar.coefs[2:4]), lag.max = 10, pacf=TRUE) #Check v[1:20] @ The matrix $\texttt{Phi.mat}$ contains the results of the Levinson-Durbin algorithm, and its diagonal elements are the PACF values match those computed by the \texttt{ARMAacf} function. Note that the PACF values are \[ \alpha_X(h) = \left\{ \begin{array}{lc} 1 & h=0 \\ \Sexpr{round(Phi.mat[1,1],4)} & h=1 \\[3pt] \Sexpr{round(Phi.mat[2,2],4)} & h=2\\[3pt] \Sexpr{round(Phi.mat[3,3],4)} & h=3 \\[3pt] \Sexpr{round(Phi.mat[4,4],4)} & h\geq 4\\ \end{array} \right. \] Also, note that the \texttt{Phi.mat} matrix has identical rows after the third. This is because for an $\AR(3)$ model, the optimal prediction is always the conditional mean. Note also that the prediction variance converges to a constant. \pagebreak \textbf{The Innovations Algorithm:} Let $\{X_t\}$ be zero-mean process with $\E[X_t^2]<\infty$ and let ${\gamma_X}(t,s)=\E[X_tX_s]$. Let \begin{align*} \widehat{X}_n=\begin{cases} 0 & \cif n=1 \\ \E[X_n|\bX_{1:(n-1)}] & \cif n \geq 2 \end{cases} \end{align*} Let $U_n=X_n-\widehat{X}_n$ be the error in prediction and let $\bs{U}_{1:n}$ be the vector of prediction errors (computed in a one-step ahead fashion). We assume that $\E[X_n|\bX_{1:(n-1)}]$ is linear in $X_1, \ldots, X_{n-1}$, and write $\bs{U}_{1:n}=\bs{A}_n\bX_{1:n}$ where $\bs{A}_n$ is given by \begin{align*} \bs{A}_n= \begin{bmatrix} 1 & 0&0 &0 &0\\[6pt] a_{1,1} & 1 & 0& 0&0\\[6pt] a_{2,2} & a_{2,1} & 1 & 0 & \vdots\\[6pt] \vdots & \vdots& \vdots& \ddots & \vdots\\[6pt] a_{n-1,n-1} &\cdots & \cdots& a_{n-1,1} &1 \end{bmatrix} \end{align*} as \begin{align*} U_1&=X_1\\U_2&=X_2-\widehat{X}_2 =X_2+a_{1,1}X_1 \\U_3&=X_3-\widehat{X}_3=X_3+a_{2,2}X_2+a_{2,1}X_1\\& \qquad \vdots \end{align*} Let $\bs{C}_n=\bs{A}_n^{-1}$, again a lower triangular matrix, where \begin{align*} \bs{C}_n=\begin{bmatrix} 1 & 0&0 &0 &0\\[6pt] \theta_{1,1} & 1 & 0& 0&0\\[6pt] \theta_{2,2} & \theta_{2,1} & 1 &0 & 0\\[6pt] \vdots & \vdots& \vdots& \ddots& \vdots\\ \theta_{n-1,n-1} &\cdots & \cdots& \theta_{n-1,1} &1 \end{bmatrix} \end{align*} Then $\bX_{1:n}=\bs{C}_n\bs{U}_{1:n}=\bs{C}_n(\bX_{1:n}-\widehat{\bX}_{1:n})$, and \[ \widehat{\bX}_{1:n} =\bX_{1:n}-\bs{U}_{1:n}=\bs{C}_n\bs{U}_{1:n}-\bs{U}_{1:n}=(\bs{C}_n-\mathbf{I} _n)\bs{U}_{1:n} = (\bs{C}_n-\mathbf{I}_n)(\bX_{1:n}-\widehat{\bX}_{1:n} )=\bs{\Theta}_n(\bX_{1:n} -\widehat {\bX }_{1:n}) \] where \begin{align*} \bs{\Theta}_n = \bs{C}_n - \mathbf{I}_n = \begin{bmatrix} 0 & 0&0 &0 &0\\[6pt] \theta_{1,1} & 0 & 0& 0&0\\[6pt] \theta_{2,2} & \theta_{2,1} & 0 &0 & 0\\[6pt] \vdots & \vdots& \vdots& \ddots& \vdots\\ \theta_{n-1,n-1} &\cdots & \cdots& \theta_{n-1,1} &0 \end{bmatrix} \end{align*} and so \begin{align*} \widehat{X}_{n+1} = \begin{cases} 0 & \cif n=1 \\[6pt] \sum\limits_{i=1}^n \theta_{n,i} (X_{n+1-i}-\widehat{X}_{n+1-i}) & \cif n \geq 2 \end{cases} \end{align*} \textit{i.e.} $\widehat{X}_{n+1}$ is a linear combination of $X_1-\widehat{X}_1$, $X_2-\widehat{X}_2$, \ldots, $X_n-\widehat{X}_n$, the one-step ahead prediction errors which are \emph{uncorrelated}. \pagebreak For the general (non-stationary) case: \begin{itemize} \item Initialize \[ v_0= {\widetilde\gamma_X}(1,1)=\E[X_1^2] \] \item Recursion at step $i=1,2,3,\ldots$: \begin{align*} \theta_{i,i-k}=\frac{1}{v_k} \left[ {\widetilde\gamma_X}(i+1,k+1)- \sum_{j=0}^{k-1} \theta_{k,k-j}\theta_{i,i-j}v_j\right] \end{align*} for $0 \leq k >= Gamma.X<-Variance*ARMAacf(ar=Re(-ar.coefs[2:4]), lag.max = n+1) v<-rep(0,n+1) v[0+1]<-Gamma.X[0+1] #Index 0 Theta.mat<-matrix(0,n+1,n+1) Theta.mat[1,1]<-Gamma.X[1+1]/v[0+1] #Index i=1 v[1+1]<-Gamma.X[0+1]-Theta.mat[1,1]^2*v[0+1] for(i in 2:n){ #Index i=2,3,... for(k in 0:(i-1)){ jvals<-0:(k-1) sumth<-sum(Theta.mat[k,k-jvals]*Theta.mat[i,i-jvals]*v[jvals+1]) Theta.mat[i,i-k] = (Gamma.X[(i-k)+1] - sumth)/v[k+1] } jvals<-0:(i-1) v[i+1]<-(Gamma.X[0+1] - sum(Theta.mat[i,i-jvals]*Theta.mat[i,i-jvals]*v[jvals+1])) } v[1:20] Theta.mat[1:5,1:5] @ For a real series, we can see how the forecast behaves: <>= set.seed(26) n<-200 X<-arima.sim(n=n+1,model=list(ar=Re(-ar.coefs[2:4]))) Xhat<-rep(0,n+1) for(ni in 1:n){ Xhat[ni+1]<-Theta.mat[ni,1:ni] %*% rev(X[1:ni] - Xhat[1:ni]) } par(mar=c(4,6,0,1)) plot(X,type='l',ylim=range(-12,12)) points(X,pch=19,cex=0.5) lines(Xhat,type='l',lty=2,col='red') points(Xhat,pch=19,cex=0.5,col='red') @ The one-step ahead forecasts seem largely to be accurate. <>= par(mar=c(4,6,0,1)) plot(as.numeric(X-Xhat),type='l',ylim=range(-3.5,3.5),ylab=expression(X-hat(X))) points(as.numeric(X-Xhat),pch=19,cex=0.5) acf(as.numeric(X-Xhat),lag.max=20,ylim=range(-1,1)) @ \end{document}