% reduce var to state-space using svd
% controls c=S*states
% states states=PP*states(-1)+Q2;

function [S,PP,Q1,Q2,S1,S2,P21,P22,N,states]=var2stsp(PHI,nsh,states)
S1=[];S=[];
PP=[];
S2=[];P21=[];P22=[];N=[];
P=ss_compact(PHI,nsh);
others=setdiff([1:size(P,1)]',states);
ns=length(states);
no=length(others);
P_s=P([states,others],[states,others]);
[U,Sv,V] = svd(P_s);
tmp=round(abs(diag(Sv))*1e10)*1e-10;
zer=find(tmp==0);
nonzer=setdiff([1:length(tmp)]',zer);
Vinv=(U');
V21=Vinv(ns+1:end,1:ns);
V22=Vinv(ns+1:end,ns+1:end);
if rcond(V22)<1e-12
    disp('non-invertible matrix')
% return
end
S=-(V22)\V21;

% s=(P_s(1:ns,1:ns)+P_s(1:ns,ns+1:end)*S)*s(-1)+M*epsilon

PP=(P_s(1:ns,1:ns)+P_s(1:ns,ns+1:end)*S);
P21=PP(nsh+1:end,1:nsh);
P22=PP(nsh+1:end,nsh+1:end);
S1=S(:,1:nsh)+S(:,nsh+1:end)*P21;
S2=S(:,nsh+1:end)*P22;
N=P_s(1:nsh,1:nsh);
M2=[eye(nsh);zeros(size(P,1)-nsh,nsh)];
Q=Vinv*M2;
Q1=Q(1:ns,:);
Q2=Q(ns+1:end,:);
return

