3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Michael A Sep 10 '15 at 21:16
  • 1
    $\begingroup$ 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. $\endgroup$ – Creosote Sep 13 '15 at 14:06
  • $\begingroup$ 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. $\endgroup$ – Creosote Sep 13 '15 at 14:18
1
+50
$\begingroup$

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. a graph

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")
add_points = function(Y) {
  points(1:t,Y[1:t],col="grey",type="l")
  points(1:t_burn,Y[1:t_burn],col="red",pch=20,type="o")
  points((t_burn+1):t,Y[(t_burn+1):t],col="blue",pch=20,type="o")
}
add_points(Y1)
add_points(Y2)
$\endgroup$
  • $\begingroup$ This is quite helpful; thank you. Your explanation of the decay makes perfect sense, and I'll look over the code over the next couple of days. $\endgroup$ – Michael A Sep 14 '15 at 13:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.