function [RHO,BETA,V] = car_simulation(rho,v,b,X,W,sims,Eig)

%SAR_SIM simulates a nearest-neighbor autoregressive process given 
%        beta, X, W, and sig = 1.
%
%INPUTS: (i)     rho  = autoregressive parameter
%        (ii)    v  = variance of errors
%        (iii)   b    = (k+1) beta vector
%        (iv)    X    = (n x k) data matrix
%        (v)     W    = (n x n) weight matrix   
%        (vi)   sims  = number of simulated data sets
%        (vii)   Eig  = (optional) eigenvalues of W 
%
%OUTPUT: Y = (n x sims)-matrix of simulated data sets.


%Augment X with unit vector

[n,k] = size(X);

XX(:,1) = ones(n,1) ;

XX(:,2:k+1) = X ;

%Preliminary Computations 

B = inv(eye(n) - rho*W) ;

mu = XX*b;

sig = sqrt(v);

A = chol(B)';

Y = zeros(n,sims);

RHO = zeros(sims,1);

V = zeros(sims,1);

BETA = zeros(k+1,sims);

%Simulate data

for i = 1:sims
   
   z = normrnd(0,1,n,1) ;
   
   y = mu + sig*A*z;
   
   [rho,b,v] = sar_sim(y,X,W,Eig);
   
   Y(:,i) = y;
   
   RHO(i) = rho;
   
   V(i) = v;
   
   BETA(:,i) = b;
   
end

rho = mean(RHO);

disp(['rho_mean = ',num2str(rho)]);

v = mean(V);

disp(['v_mean = ',num2str(v)]);


b = mean(BETA')