function [v]=arsim(w,A,C,n,ndisc)
%ARSIM	Simulation of AR process.	
%
%	v=ARSIM(w,A,C,n) simulates n time steps of the AR(p)
%	process
%
%        v(k,:)' = w' + A1*v(k-1,:)' +...+ Ap*v(k-p,:)' + eta(k,:)', 
%
%	where A=[A1 ... Ap] is the coefficient matrix and w is a 
% 	vector of intercept terms that is included to allow for a
%	nonzero mean of the process. The vectors eta(k,:) are
%	independent Gaussian noise vectors with zero mean and
%	covariance matrix C. 
%
%	The simulation requires p vectors of initial values, all of
%	which are taken to equal the mean value of the process. (The
%	process mean is calculated from the parameters A and w.) To
%	avoid possible spin-up effects, the first 10^3 time steps are
%	discarded. Alternatively, ARSIM(w,A,C,n,ndisc) discards the
%	first ndisc time steps. 

%	Modified 14-Feb-97
%	Author: Tapio Schneider
%		AOS Program, Princeton University
%		Princeton, NJ 08544-0710
%		tapio@splash.Princeton.EDU


m = size(C,1);  			% dimension of state vectors 
p = size(A,2)/m;  			% order of process

if (p ~= round(p)) 
   error('Bad arguments.'); 
end

if (length(w) ~= m | min(size(w)) ~= 1)
   error('Dimensions of arguments are incompatible with each other.')
end 

w = w(:)';  				% force w to be row vector

% Discard the first ndisc time steps; if ndisc is not given as input
% argument, use default
if (nargin < 5) 
   ndisc = 10^3; 
end

% Compute Cholesky factor of covariance matrix C
R = chol(C);  				% R upper triangular

% Get ndisc+n independent Gaussian pseudo-random vectors with 
% covariance matrix R'*R
randvec = randn([ndisc+n,m])*R;

% Add intercept vector to random vectors
randvec = randvec+ones(ndisc+n,1)*w;

% Get transpose of system matrix A (use transpose in simulation because 
% we want to obtain the states as row vectors)
AT = A';

% Initialize state vectors
u = zeros(p+ndisc+n,m);

% Take the p initial values of the simulation to equal the process mean, 
% which is calculated from the parameters A and w

if any(w)
%  Process has nonzero mean    mval = inv(B)*w'    where 
%             B = eye(m) - A1 -... - Ap; 
%  Assemble B
   B = eye(m);
   for j=1:p
       B = B - A(:, (j-1)*m+1:j*m);
   end
%  Get mean value of process
   mval = w / B';
   x = ones(p,1)*mval;
else
%  Process has zero mean
   x = zeros(p,m); 
end

% The vector x serves as intermediate storage of the optimal forecast 
% of the next state given the p previous ones and is now initialized with 
% the process mean.

% Simulate n+ndisc observations
for k=p+1:ndisc+n; 
    for j=1:p;
        x(j,:) = u(k-j,:)*AT((j-1)*m+1:j*m,:);
    end
    u(k,:) = sum(x)+randvec(k,:);
end
v = u(ndisc+1:ndisc+n,:); % return only the last n state vectors





