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

A: 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:


*

*Generate eps, a vector of 1100 draws from the Gaussian distribution, $\epsilon_t$;

*Initialize the first 3 elements in the vector x with zeros: x[1] = x[2] = x[3] = 0.

*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.
A: 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.
A: Use the rGARMA function in the ts.extend package
If you don't mind using R instead of MATLAB, you can generate random vectors from any stationary Gaussian ARMA model (including AR models) using the ts.extend package.  This package generates random vectors directly form the multivariate normal distribution using the computed autocorrelation matrix for the random vector, so it gives random vectors from the exact distribution and does not require "burn-in" iterations.  Here is an example from an $\text{AR}(3)$ model.
#Load the package
library(ts.extend)

#Set parameters
AR <- c(0.8,  -0.3, 0.1)
m  <- 50

#Generate n = 12 random vectors from this model
set.seed(1)
SERIES <- rGARMA(n = 12, m = m, ar = AR)

#Plot the series using ggplot2 graphics
library(ggplot2)
plot(SERIES)


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