function [shortf,invDf,fhat] = simstate_chanjel(yt,lam,phi,psi,sigeps,sigeta,psilag,philag,nf)

nvar=size(yt,2);
nobs=size(yt,1);

%settings for chan&jeliazkov
T=nobs;
Tn=nobs*nvar;
Tf=(nobs+max(philag,psilag));
Tq=(nobs+max(philag,psilag))*nf;

% enter system matrices for state space model 
% yt = hmat*zt + e 
% zt = ct + fmat*zt-1 + u 

%construct the banded lambda matrix
hmat = lam;
if psilag
for iind=1:psilag;
   hmat = [-lam.*kron(psi(:,iind),ones(1,nf)) hmat];
end
end
% if size(hmat,2)<zdim
%     hmat=[zeros(nvar,zdim-size(hmat,2)) hmat];
% end
hmat=repmat(hmat,nobs,1)';

ind_r=kron(ones(1,nvar),1:nf*(psilag+1));
ind_c=kron(1:nvar*nobs,ones(1,nf*(psilag+1)));
ind_r=[ind_r(ones(1,nobs),:)+kron([0:nobs-1]'*nf,ones(1,nvar*nf*(psilag+1)))]';
ind_r=ind_r(:)';
hmat_sp=sparse(ind_r,ind_c,hmat(:),(nobs+psilag)*nf,nvar*nobs)';
if philag>psilag
    hmat_sp=[sparse(nvar*nobs,nf*(philag-psilag)) hmat_sp];
end

%construct the banded phi matrix
ind_f=kron(ones(1,philag),(nf*(philag-1)+1):nf*philag)-kron([0:philag-1]*nf,ones(1,nf));
hmat=[-phi(:,ind_f) eye(nf)];
hmat=repmat(hmat,Tf-philag,1)';

ind_r=kron(ones(1,nf),1:nf*(philag+1));
ind_c=kron(1:nf*(Tf-philag),ones(1,nf*(philag+1)));
ind_r=[ind_r(ones(1,Tf-philag),:)+kron([0:Tf-philag-1]'*nf,ones(1,nf^2*(philag+1)))]';
ind_r=ind_r(:)';
fmat_sp=[speye(Tq,philag*nf) sparse(ind_r,ind_c,hmat(:),Tq,Tq-philag*nf)]';

%initial state variance and factor innovation variance matrix
sigeta0=speye((nf*philag)^2);
if philag>0
    Sig=[sigeta,sparse(nf,nf*(philag-1));sparse(nf*(philag-1),nf*philag)];
    sigeta0=sigeta0-kron([phi;speye(nf*(philag-1),nf*philag)],[phi;speye(nf*(philag-1),nf*philag)]);
    sigeta0=reshape(inv(sigeta0)*Sig(:),nf*philag,nf*philag);
    sigeta0=speye(nf,nf*philag)*sigeta0*speye(nf*philag,nf);

    Sinv0=inv(sigeta0);
    Sinv=inv(sigeta);
    Sinv_sp=sparse([kron(speye(philag),Sinv0) sparse(philag*nf,(Tf-philag)*nf);...
        sparse((Tf-philag)*nf,philag*nf) kron(speye(Tf-philag),Sinv)]);
else
    Sinv=inv(sigeta);
    Sinv_sp=sparse([kron(speye(Tf-philag),Sinv)]);
end

F0=fmat_sp'*Sinv_sp*fmat_sp; %prior precision

postinv=hmat_sp'*kron(speye(nobs),spdiags(1./sigeps,0,nvar,nvar));
invDf = F0 + postinv*hmat_sp;  % posterior precision
df = postinv*reshape(yt',Tn,1);
Cf = chol(invDf);
fhat = Cf\(Cf'\df);          % posterior mean
f = fhat + Cf\randn(Tq,1);   
shortf = reshape(f,nf,Tf)';

   
    