# How to simulate a third order AR model

I'm trying to understand AR models but it's getting pretty difficult for me. Just wanted to ask you some hints on how to simulate an AR(3) model driven by a zero mean WN for 1000 values in Matlab, without using any built function. A recommendation on a good source for understanding this would work as well.

• Both given answers are excellent, but note that there are restrictions on coefficients on AR(3) model, i.e. not all parameter combinations give valid AR(3) model. – mpiktas Apr 13 '15 at 10:35

Its not common to give away the answer for a question but I had this code ready so I can just as well post it together with an explanation of what each line does. This code does however not discard the first 100 observations so that is something you can do on your own.

%Simulate AR(3)
T     = 1000;                         %Set how many observations you need
y     = ones(T,1);                    %Create a vector of dim Tx1 to store the simulations in
y(1)  = 1;                            %Set the first obs. to 1
y(2)  = 0.5;                          %Set the second obs. to 0.5
y(3)  = 1.5;                          %Set the third obs. to 1.5
rho1  = 0.2;                          %Set the value of rho1 (coefficient on y(t-1))
rho2  = 0.2;                          %Set the value of rho2 (coefficient on y(t-2))
rho3  = 0.1;                          %Set the value of rho3 (coefficient on y(t-3))
sigma = 1;                            %Set the value of the s.d. of the error term
mu_e  = 0;                            %Set the value of the mean of the error term
eps   = normrnd(mu_e, sigma, T, 1);   %Creat a vector of normal random numbers with mean, mu_e and s.d. sigma. Dimension is Tx1

for t=4:1000;                         %Start the loop running from obs. 4 to 1000
y(t) = rho1*y(t-1) + rho2*y(t-2) + rho2*y(t-3) + eps(t);    %The AR(3) model
end

%Plot the series
figure
plot(y);
title('AR(3)');
xlabel('t')
ylabel('y(t)')


Let's define the third order autoregressive model, AR(3), as follows:

$$x_t = \alpha_1 x_{t-1} + \alpha_2 x_{t-2} + \alpha_3 x_{t-3} + \epsilon_t \,,\;\; \epsilon_t \sim NID(0, \sigma^2) \,, \hbox{ for } t=1,2,\dots, n\,.$$

A straightforward way to generate data from the equation above is by means of a loop. These are the steps and some pseudocode:

1. Generate eps, a vector of 1100 draws from the Gaussian distribution, $\epsilon_t$;
2. Initialize the first 3 elements in the vector x with zeros: x[1] = x[2] = x[3] = 0.
3. Run the loop: for i=4 to i = 1100 do x[i] = a1 * x[i-1] + a2 * x[i-2] + a3 * x[i-3] + eps[i] and discard the first 100 observations in x to reduce the effect of initialization in step 2.

The above is a do-it-yourself solution that may be interesting for pedagogical purposes. For production, it is better to use specialized functions already available in software packages. For example, in Matlab see this piece of documentation.

Write the process as

$$\begin{pmatrix} y_t \\ y_{t-1} \\ y_{t-2} \end{pmatrix} = \underbrace{\begin{pmatrix}\phi_1 & \phi_2 & \phi_3 \\ 1 & 0 & 0 \\ 0 & 1 & 0\end{pmatrix}}_F \begin{pmatrix} y_{t-1} \\ y_{t-2} \\ y_{t-3} \end{pmatrix} + \begin{pmatrix} \epsilon_t \\ 0 \\ 0 \end{pmatrix}, \qquad \epsilon_t\sim N(0,\sigma^2)$$

Then save some computations by finding the time invariante covariance matrix which satisfies $$Q_0 = \begin{pmatrix}\sigma^2 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} + FQ_0F^\top$$

assuming that the process is stationary. Use this to sample the first index. R code to do this is something like

# define loadings
p <- 3L
F. <- matrix(0., p, p)
F.[1, ] <- c(0.1, 0.5, 0.25)                      # phi_1, phi_2, phi_3
F.[2:p, 1:(p - 1L)] <- diag(p - 1L)

# get time-invariant covariance matrix
sig <- .2                                         # sigma
Q_0 <- local({
Q <- matrix(c(sig^2, numeric(p * p - 1L)), p, p)
eg <- eigen(F.)
las <- eg$$values U <- eg$$vectors
T. <- solve(U, t(solve(U, Q)))
Z <- T./(1 - tcrossprod(las))
Re(tcrossprod(U %*% Z, U))
})

# confirmation that Q_0 is correct
Q <- matrix(c(sig^2, numeric(p * p - 1L)), p, p)
X <- matrix(c(sig^2, numeric(p * p - 1L)), p, p)
for(t. in 1:1000)
X <- Q + tcrossprod(F. %*% X, F.)
all.equal(Q_0, X)
#R [1] TRUE

# use Reduce to simulate and confirm estiamtes with ar
set.seed(24456488)
n_sim <- 100L
a_0 <- drop(crossprod(chol(Q_0), rnorm(p)))
sims <- Reduce(
function(x, y){
o <- F. %*% x
o[1L] <- o[1L] + y
o
}, rnorm(n_sim, sd = sig), init = a_0,
accumulate = TRUE)
sims <- sapply(sims, "[[", 1L)
ar(sims, aic = FALSE, order.max = p)
#R
#R Call:
#R ar(x = sims, aic = FALSE, order.max = p)
#R
#R Coefficients:
#R      1       2       3
#R 0.1434  0.4591  0.1112
#R
#R Order selected 3  sigma^2 estimated as  0.0452


I have not worked much with Matlab but I gather one can easily convert the above to Matlab code.

Depends on whether you accept other built-in functions which are not specifically for AR processes, but help a lot with generating one. So, firstly use Matlab's randn function to get a vector of normally distributed i.i.d. random values (with zero mean and unit variance):

e = randn(1000, 1);


Then filter this signal with an all-pole filter to get the desired AR process:

x = filter(1, [1 a], e);


Here the first 1 corresponds to the filter having no zeros, or in other words not having the moving average part (see ARMA process for more info). [1 a] is the vector of filter coefficients specifying the poles (for order of 3, we would have e.g. a = [1 2 3];).