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.
 A: 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")
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)

