AR(1) with known initial and terminal condition: how to draw the innovations?

Question

Suppose I have the following stationary $AR(1)$ process:

$$ y_{t}=\alpha_{0}+\alpha_{1}y_{t-1} + u_{t} $$

where $u_{t} \sim \mathbb{N}(0,\sigma^{2})$, with $\sigma^{2}$ known. Suppose I have an initial condition $y_{1}$ and terminal condition $y_{T}$ and I would like to simulate my process for the periods in the interim, i.e $t = 2,\dots,T-1$. After substituting recursively, I can write the final term $y_{T}$ as

$$ y_{T} = \alpha_{0} \left(\frac{1-\alpha_{1}^{T-1}}{1-\alpha_{1}} \right) + \alpha_{1}^{T-1}y_{1} + \sum_{t=0}^{T-2} \alpha_{1}^{t}u_{T-t} $$

Since $y_{T}$ and $y_{1}$ are known, this suggests that the sum of the innovations should equal

$$ \sum_{t=0}^{T-2} \alpha_{1}^{t}u_{T-t} = y_{T} - \alpha_{0} \left(\frac{1-\alpha_{1}^{T-1}}{1-\alpha_{1}} \right) - \alpha_{1}^{T-1}y_{1} $$

Now, my question is what is the correct way to draw the sequence of $u_{t}$ knowing that they are $iid$ and unconditionally distributed as a $\mathbb{N}(0,\sigma^{2})$.

The way I am trying to approach this problem is following the comment of @soakley on this post Sample random variables conditional on their sum . However, I am not sure about the proper iterative procedure.

Clearly, I could treat the $\alpha_{1}^{t}$ as weights. But I am still confused. The way I proceed is as follows:

0) Define $z_{2}=y_{T} - \alpha_{0} \left(\frac{1-\alpha_{1}^{T-1}}{1-\alpha_{1}} \right) - \alpha_{1}^{T-1}y_{1}$

1) Draw $u_{t}\sim \mathbb{N}\left( \frac{z_{t}}{T-t+1},\frac{T-t}{T-t+1} \sigma^{2} \right)$ for $t \in [2,T-1]$

2) In each step update $z_{t}$ as follows

$$ z_{t+1} = z_{t} - \alpha^{T-t} u_{t} $$

3) Let $u_{T} = z_{T}$

I would really like to know if this procedure is correct.

Answer 1

Since the overarching goal here is to generate the $\text{AR}(1)$ series values, I will show you how you can generate these directly, rather than via the error terms. (The method of generating the error terms is very similar.) The way to do this is to write out the multivariate normal distribution for the process and then use the rules for the conditional distribution for a multivariate normal. I will use more standard notation by using the form:

$$Y_t = \mu + \phi (Y_{t-1} - \mu) + u_t \quad \quad \quad u_t \sim \text{IID N}(0, \sigma^2).$$

The stationary form of this model is for vectors of $y_i$ values to be distributed by a multivariate normal distribution, with variance matrix determined by the auto-covariance function of the process. Thus, when you want conditional distributions, you can use standard results for conditioning in the multivariate normal distribution.

---

Conditional distribution given endpoints: The auto-correlation function for the $\text{AR}(1)$ process above is:

$$\rho(t) = \frac{\phi^t}{1-\phi^2} \quad \quad \quad \text{for } t \in \mathbb{Z}.$$

Thus, the marginal distribution of the vector $\mathbf{Y} = (Y_0,...,Y_T)$ is:

$$\mathbf{Y} \sim \text{N} \Bigg( \mu \mathbf{1}, \sigma^2 \mathbf{\Sigma}(\phi) \Bigg) \quad \quad \quad \mathbf{\Sigma}(\phi) \equiv \frac{1}{1-\phi^2} \begin{bmatrix} 1 & \phi & \phi^2 & \cdots & \phi^T \\ \phi & 1 & \phi & \cdots & \phi^{T-1} \\ \phi^2 & \phi & 1 & \cdots & \phi^{T-2} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \phi^T & \phi^{T-1} & \phi^{T-2} & \cdots & 1 \\ \end{bmatrix}.$$

To obtain the conditional distribution of interest, we decompose the random vector $\mathbf{Y}$ as:

$$\mathbf{y} = \begin{bmatrix} \mathbf{Y}_\text{INT} \\ \mathbf{Y}_\text{END} \end{bmatrix},$$

where $\mathbf{Y}_\text{INT} = (Y_1,...,Y_{T-1})$ and $\mathbf{Y}_\text{END} = (Y_0,Y_T)$. We likewise decompose the mean vector and variance matrix as:

$$\mu \cdot \mathbf{1} = \begin{bmatrix} \mu \cdot \mathbf{1} \\ \mu \cdot \mathbf{1} \end{bmatrix} \quad \quad \quad \mathbf{\Sigma}(\phi) = \begin{bmatrix} \mathbf{\Sigma}_\text{INT}(\phi) & \mathbf{\Sigma}_\text{CROSS}(\phi) \\ \mathbf{\Sigma}_\text{CROSS}(\phi)^\text{T} & \mathbf{\Sigma}_\text{END}(\phi) \end{bmatrix},$$

where the mean vectors are of appropriate lengths and the variance parts are:

$$\begin{aligned} \mathbf{\Sigma}_\text{INT}(\phi) &\equiv \frac{1}{1-\phi^2} \begin{bmatrix} 1 & \phi & \phi^2 & \cdots & \phi^{T-2} \\ \phi & 1 & \phi & \cdots & \phi^{T-3} \\ \phi^2 & \phi & 1 & \cdots & \phi^{T-4} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \phi^{T-2} & \phi^{T-3} & \phi^{T-4} & \cdots & 1 \\ \end{bmatrix}, \\[30pt] \mathbf{\Sigma}_\text{CROSS}(\phi) &\equiv \frac{1}{1-\phi^2} \begin{bmatrix} \phi & \phi^{T-1} \\ \phi^2 & \phi^{T-2} \\ \phi^3 & \phi^{T-3} \\ \vdots & \vdots \\ \phi^{T-1} & \phi \\ \end{bmatrix}, \\[30pt] \mathbf{\Sigma}_\text{END}(\phi) &\equiv \frac{1}{1-\phi^2} \begin{bmatrix} 1 & \phi^{T} \\ \phi^{T} & 1 \\ \end{bmatrix}. \\[6pt] \end{aligned}$$

Now, using the conditional density rule for the multivariate normal distribution, we have:

$$p(\mathbf{y}_\text{INT}| \mathbf{y}_\text{END}) = \text{N}(\mathbf{y}_\text{INT}| \boldsymbol{\mu}_*, \boldsymbol{\Sigma}_*),$$

where:

$$\begin{equation} \begin{aligned} \boldsymbol{\mu}_* &\equiv \mu \cdot \mathbf{1} + \boldsymbol{\Sigma}_\text{CROSS}(\phi) \boldsymbol{\Sigma}_\text{END}(\phi)^{-1} (\boldsymbol{y}_\text{END} - \mu \cdot \mathbf{1}), \\[8pt] \boldsymbol{\Sigma}_* &\equiv \boldsymbol{\Sigma}_\text{INT}(\phi) - \boldsymbol{\Sigma}_\text{CROSS}(\phi) \boldsymbol{\Sigma}_\text{END}(\phi)^{-1} \boldsymbol{\Sigma}_\text{CROSS}(\phi)^\text{T}. \\[6pt] \end{aligned} \end{equation}$$

This gives you the conditional distribution of the observable values in the process. If you particularly want the error terms, you can extract them from the observed values using the recursive equation for the process.

Answer 2

Here is my revised and (hopefully) final attempt.

Suppose we have the AR(1) process: $y_{t}= \alpha \times y_{t-1} + \epsilon_t$ and we know $y_1$ and $y_{T}$ and want to figure out the innovations that will start at $t = 2$ and bring us to $y_{T}$.

So, suppose that there is an innovation at time $t = 2$ and it is $\epsilon_{2}$ and similarly a new $\epsilon_{3}$ at $t= 3$ and the same thing upto and including $t = {T-1}$.

Then, what will this result in as far as $y_{2}, y_{3}, \ldots y_{T-1}$.

Well, then we have ( don't ask what I was doing in my earlier attempts ):

$y_{2} = \alpha \times y_{1} + \epsilon_{2} $

$y_{3} = \alpha^2 y_{1} + \alpha \times \epsilon_{2} + \epsilon_{3} $

$y_{4} = \alpha^3 y_{1} + \alpha^2 \times \epsilon_{2} + \alpha \epsilon_{3} + \epsilon_{4} $

So, at time $t={T-1}$, we have

(1) $y_{(T-1)} = \alpha^{T-1} y_{T} + \sum_{i=2}^{T-1} \epsilon_{i} \alpha^{T-i}$

But, since $y_{T} = \alpha \times y_{T-1} + \epsilon_{T}$, this leads to

(3) $ y_{T} = \alpha^{T} y_{1} + \alpha \sum_{i=2}^{T-1} \epsilon_{i} \alpha^{T-i} + \epsilon_{T}$

Notice that the last term is obtained by taking the term at time t=T out of the summation. This makes it very easy to solve for $\epsilon_{T}$ algebraically which leads to:

(4) $\epsilon_{T} = y_{T} - \alpha^{T} y_{1} - \alpha \sum_{i=2}^{T-1} \epsilon_{i} \alpha^{T-i}$.

Everything on the RHS of (4) is known which means that $\epsilon_{T}$ is a function of variables whose values are all known at time $t = T$.

And my sincere apologies for all the noise and mistakes. Your expressions were totally correct so maybe the procedure that you describe and the one I desribe here are equivalent ? I would think that the values of $\epsilon_{t}$ from $t=2$ to $t=T-1$ can be anything as long as $\epsilon_{T}$ is the deterministic function of them given by (4). But I've lost all confidence in my ability to solve this problem correctly so correct me if you're thinking is different. :).

Answer 3

Use the rGARMA function in the ts.extend package

(Note: This answer augments my other theory-based answer, the method for which has now been programmed into the present package.)

You can generate random vectors from the marginal or conditional distributions of any stationary Gaussian ARMA model using the ts.extend package. This package generates random vectors directly form the multivariate normal distribution using the computed autocorrelation matrix for the random vector, so it gives random vectors from the exact distribution and does not require "burn-in" iterations. It also allows you to impose conditioning values and in this case the package generates from the resulting conditional distribution.

Here is an example from an $\text{AR}(1)$ model where you generate values at the terminal points at either end of the process. (Here I will use a model with parameters $\alpha_0 = 0$, $\alpha_1 = 0.4$ and $\sigma=1$. I will generate twelve random vectors from this model each with $m=30$ points with terminal conditioning points $y_1=-3$ and $y_m=3$.)

#Load the package
library(ts.extend)

#Set parameters
AR <- 0.4
m  <- 30
y1 <- -3
ym <- 3

#Set the conditioning vector
CONDVALS    <- rep(NA, m)
CONDVALS[1] <- y1
CONDVALS[m] <- ym

#Generate n = 12 random vectors from this model
set.seed(1)
SERIES <- rGARMA(n = 12, m = m, ar = AR, condvals = CONDVALS)

#Plot the series using ggplot2 graphics
library(ggplot2)
plot(SERIES)

As you can see from the image below, each of the plots shows a time-series vector that has the stipulated starting and ending points. The interior points are generated from the conditional distribution for the process. As you can see in the code, the syntax for this is quite simple --- all you need to do is to create the vector CONDVALS with the stipulated conditioning values and add this as an argument for the function.