% SIM_DAT - Simulate date from the model solution
%
% usage
%
% [xsim, ysim, shock] = sim_dat(gx,hx,eta,nper,ndrop)

function [xsim, ysim, shock] = sim_dat(gx,hx,eta,nper,ndrop, x0)
ny = size(gx,1);
nx = size(hx,1);
neps = size(eta,2);

shock = randn(neps,nper);
ysim = zeros(ny,nper);
xsim = zeros(nx,nper);

%Initial state
if nargin>5
    xsim(:,1) = x0;
end

for j = 1:nper-1
    ysim(:,j) = gx*xsim(:,j);
    xsim(:,j+1) = hx*xsim(:,j)+eta*shock(:,j+1);
end
ysim(:,j+1) = gx*xsim(:,j+1);


%Drop burn-in periods
xsim = xsim(:,ndrop+1:end);
ysim = ysim(:,ndrop+1:end);
shock = shock(:,ndrop+1:end);
