# Simulate from an AR(1) process where the parameters vary over time

I would like to simulate from an AR(1) process, where the $\rho$ parameter in the process of the form:

$X_t = \rho X_{t-1} + \epsilon_t$

varies over time. The path of $\rho$ and $var(\epsilon_t)$ over time is pre-determined by some other process; for argument's sake here, let us say that the parameters are given beforehand in the lists:

listRho = RandomReal[1, {100}];
listSigma = RandomReal[1, {100}];


I know how to simulate from an AR(1) process (for example where $\rho=0.5$ and $var(\epsilon_t) = 1$) using:

 RandomFunction[ARProcess[0, {0.5}, 1]


However, does anyone know whether I can use the in-built functions, or otherwise, to allow $\rho$ and $var(\epsilon_t)$ to vary at each time step?

Update: I can simulate from an AR(1) manually using FoldList:

FoldList[0.5 #1 + #2 &, 0, RandomVariate[NormalDistribution[0, 1], {100}]]


I have tried to allow the parameters of this process to vary over time using:

Reap[t = 1; FoldList[t++; Sow[t]; listRho[[t]] #1 + #2 &, 0, RandomVariate[NormalDistribution[0, listSigma[[t]]], {100}]]]


Where I use Reap and Sow to try to check whether $t$ is incrementing, but it does not appear to be.

Anyone know how I can do this?

Best,

Ben

First let's define the needed variables:

n = 100;
ρ = RandomReal[1, {n}];
σ = RandomReal[1, {n}];
noise = RandomVariate[NormalDistribution[], n]*σ;


Note that in order to get the noise, I multiply a list of noises drawn from a standard normal by the standard deviations. This is much faster than calling RandomVariate hundreds of times. Now here's the actual data generation:

data = Rest@FoldList[{#1, 1}.#2 &, 0, Transpose[{ρ, noise}]]


First we transpose the list of parameters and noise values together. This allows us to pass each pair of $(\rho_i,\epsilon_i)$ to the function in FoldList, without resorting to indexing. We initialize with 0. The function itself takes the dot product {#1, 1}.#2, which is simply {x[i-1], 1}.{ρ[i],ε[i]}; so we get x[i] = x[i-1]ρ[i] + ε[i]. If you want to add a constant (time-varying or not), simply add it to the noise vector.

• Thanks for that. I really appreciate the time taken to write the answer! Much better than mine. Best, Ben – ben18785 Feb 25 '15 at 18:22
• I did a quick timing comparison when I saw your answer (I started writing this answer before yours got posted). Mine is around 4.5 times faster generating 1M datapoints. – 2012rcampion Feb 25 '15 at 18:27
• Nice. Thanks for doing that. I concede! Best, Ben – ben18785 Feb 25 '15 at 18:28

I have figured out how to do this, and have checked that all is working as should be, using Reap and Sow.

First, we generate the lists, and use the $\sigma$ one to generate the error list:

listSigma = RandomReal[1, {10}];
listRho = RandomReal[1, {10}];
lError = RandomVariate[NormalDistribution[0, #]] & /@ listSigma;


I have included the last bit with Thread, so that we can compare it with the entries from Reaping and Sowing, when we run the AR(1) process below:

FNew[rho_, xPast_, error_] := Module[{}, Sow[{rho, error}]; t++; rho xPast + error]
Reap[t = 1; FoldList[FNew[listRho[[t]], #1, #2] &, 0, lError]][[2]]
Reap[t = 1; FoldList[FNew[listRho[[t]], #1, #2] &, 0, lError]][[1]]


The first reaping/sowing yields the same parameters as from the threading, indicating that the parameters are being updated at each time step. The latter reap/sow yields the AR(1) process.

Hope this helps someone else!

Best,

Ben