# How do I quantify the decay in the initial condition of an AR process?

I'm working on basic code that generates data from AR and VAR processes; the code generates enough observations to dampen the effect of the initial conditions. For example, if I want to generate 30 observations, the code generates 300 observations and only returns the final 30.

In an answer to my question about this code , a user recommended that instead of generating $$N$$ observations and only returning the last $$T$$, I pass a parameter that specifies

the desired decay of the initial condition.

I'm not sure how to interpret this. Take a simple univariate AR(1) process as an example:

\begin{align} y_t = \rho y_{t-1} + u_t \end{align} where $$\lvert \rho \rvert < 1$$, $$u \sim N(0, \sigma_u^2)$$, and $$y_0 = 0$$. What does it mean to talk about the desired decay of this initial condition? I couldn't find any references that discussed this. For example, is this referring to calculating some number of observations after which the initial condition has decayed by, say, 90%? What exactly does that mean?

The actual processes I'm working with are slightly more involved, e.g.

\begin{align} \begin{bmatrix} y_{1,t} \\ y_{2,t} \end{bmatrix} &= \begin{bmatrix} 0.02 \\ 0.03 \end{bmatrix} + \begin{bmatrix} 0.5 & 0.1 \\ 0.4 & 0.5 \end{bmatrix} \begin{bmatrix} y_{1,t-1} \\ y_{2,t-1} \end{bmatrix} + \begin{bmatrix} 0 & 0 \\ 0.25 & 0 \end{bmatrix} \begin{bmatrix} y_{1,t-2} \\ y_{2,t-2} \end{bmatrix} + \begin{bmatrix} u_{1,t} \\ u_{2,t} \end{bmatrix} \end{align}

but I'm just trying to get a general idea for now.

• I'm starting to wonder if there are any references about this concept, if it's even a valid concept. After weeks of searching and no responses either on here or on Code Review, I'm increasingly skeptical. – Michael A Sep 10 '15 at 21:16
• I think the user is just commenting on the choice of $N=10T$ in your code. There, "10" is a magic constant. I think he's saying you should instead calculate $N$ as a function of the desired "decay" of the initial condition. For example, in your AR(1) process, $y_t = \rho^{t-1}y_0 + f(\mathbf{u}_{1,\ldots,t})$. Waiting until $t=1+\log(0.01)/\log(\rho)$ would give you $y_t=0.01y_0+f(\mathbf{u}_{1,\ldots,t})$, washing out the effect of $y_0$ to one hundreth of its value, for example. – Creosote Sep 13 '15 at 14:06
• For more complicated processes, rather than solve analytically, you might want to run a second process with no noise term in the recurrence relation. Keep stepping both processes until the variables in the no-noise process appear to have converged to steady-state values. – Creosote Sep 13 '15 at 14:18

I just remembered I'm supposed to give answers as an Answer, not a Comment, so I'll use this opportunity to offer some code in R that simulates your example process with a dynamically-determined burn-in period. Red points are burn-in; blue are steady-state.

In case you can't run R, here's an example of the output.

And if you can run R, here's the code.

# parameters of the process
# rho is 0.99; each u_t is N(0,0.02**2)
rho = 0.99
u = function() rnorm(1,mean=0,sd=0.02)

# Y is the actual process, Z the no-noise cousin
# start each process at 1,1
Y1=Y2=Z1=Z2=c(1,1)
t=2

# say a sequence is "converged" if the last three
# terms are within a factor of 1+-1e-6 of each other
converged = function(Z) {
n = length(Z)
if (n<3) F else { L=log(Z[(n-2):n]); mu=mean(L); max(abs(L-mu))<1e-6; }
}

# step until converged
while (!converged(Z1) | !converged(Z2)) {
t = t + 1
Y1[t] = 0.02 + 0.5*Y1[t-1] + 0.1*Y2[t-1] + u()
Z1[t] = 0.02 + 0.5*Z1[t-1] + 0.1*Z2[t-1] # no noise
Y2[t] = 0.03 + 0.4*Y1[t-1] + 0.5*Y2[t-1] + 0.25*Y1[t-2] + u()
Z2[t] = 0.03 + 0.4*Z1[t-1] + 0.5*Z2[t-1] + 0.25*Z1[t-2] # no noise
}

# step the same amount again for "steady state" data
t_burn = t
while (t <= 2*t_burn) {
t = t + 1
Y1[t] = 0.02 + 0.5*Y1[t-1] + 0.1*Y2[t-1] + u()
Y2[t] = 0.03 + 0.4*Y1[t-1] + 0.5*Y2[t-1] + 0.25*Y1[t-2] + u()
}

# plot the results; red=burn-in period, blue="steady state" data
plot(c(1,t),range(c(Y1,Y2)),type="n",xlab="t",ylab="Y")