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?




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.

| improve this answer | |
  • $\begingroup$ Thanks for that. I really appreciate the time taken to write the answer! Much better than mine. Best, Ben $\endgroup$ – ben18785 Feb 25 '15 at 18:22
  • 1
    $\begingroup$ 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. $\endgroup$ – 2012rcampion Feb 25 '15 at 18:27
  • $\begingroup$ Nice. Thanks for doing that. I concede! Best, Ben $\endgroup$ – 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;
Thread[{listRho, lError}]

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!



| improve this answer | |

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.