clear all
close all
randn('seed',1999)

irf_length=20;
namemod=input('Give name of model without extension  ==>> ','s');%'neoc_grwt_dynare_exp';
eval(['load matrices_',namemod]);
M_=struct('fname',namemod');
M_.endo_names=matrices.varyt;
M_.param_names=matrices.parnames;
eval(['y_ss=',namemod,'_steadystate;'])
y_ss=y_ss(end:-1:1);
namepar=input('Give name of parameter file without extension  ==>> ','s');
eval(namepar);%param_neoc_grwt_dynare;
nsh=size(matrices.shocky,1);
neq=size(matrices.varyt,1);

disp('----------------------------------------------------------')
disp('List of shocks')
disp([num2str([1:nsh]'),matrices.shocky])%(matrices.sortinnov,:)
CSigma=0;
while size(CSigma,1)~=nsh && size(CSigma,2)~=nsh
    fprintf('Need %i dimensional vector or matrix\n',nsh)
    CSigma=input('Give the SQRT(covariance matrix) of the innovations to be simulated ==> ');
end
disp('----------------------------------------------------------')
if any(size(CSigma)==1)
    CSigma=diag(CSigma);
end
Sigma_mx=CSigma^2;

do_irf=input('Compute IRFs? (1/0[default])');
if do_irf
    inno=input('Give # of shock for IRF ==> ');
end


% evaluate steady state variables
for jj=1:neq
    eval(['global ',regexprep(M_.endo_names(jj,:),'\s*',''),'_ss;']);
end

neworder=[nsh:-1:1,nsh+1:neq];

eval([namemod,'_A1']);
eval([namemod,'_A2']);
eval([namemod,'_A3']);
eval([namemod,'_X']);
eval([namemod,'_M']);

% Reorder the model to have shocks first
A1=A1(end:-1:1,end:-1:1);
A2=A2(end:-1:1,end:-1:1);
A3=A3(end:-1:1,end:-1:1);
%

%% reordering the order of the shocks to be consistent with RE solution
A1=A1(neworder,neworder);
A2=A2(neworder,neworder);
A3=A3(neworder,neworder);

%%
if ~exist('Sigma_mx','var')
    stand_dev=1/100;
    
    Sigma_mx=diag(ones(size(Xmx,1),1))*stand_dev;
end

% solution is y=[Pi11 Pi21][x;s]; s=[Xmx 0;PHI1,PHI2];
% [Pi11,Pi21,PHI1,PHI2]=solver_muc(namemod, matrices);
[Pi,PHI,AA,BB,numpre,aqz,bqz,Qqz,Zqz,inz22,algo]=firstorder_stuff(namemod,matrices,0);
[vv,covy]=variance(PHI,Pi,[Msh;zeros(neq,nsh)],Sigma_mx,3000,0,0);
%% Purge constant states
const_st=find(abs(vv)<1e-10);
C=zeros(size(AA,1),nsh);%set dummy shocks to 2 as 1 could generate problems
C(neq-nsh+1:neq,1:nsh)=eye(nsh);%
[TT,Tm,war]=testsol(PHI,Pi,nsh,size(AA,1),AA,BB,C, Xmx);

new_varnames=[M_.endo_names(nsh:-1:1,:);M_.endo_names(nsh+1:end,:)];
nstat=neq;
njump=find(sum(abs(A1),1)>1e-10);
PHI2=PHI;%([1:nsh,2*nsh+1:end],[1:nsh,2*nsh+1:end]);

Pi2=Pi;%(:,[1:nsh,2*nsh+1:end]);%res2=max(max(abs(A1*RET^2+A2*RET+A3)));
%% SHIFTING OF STATES TO GET ALL ENDOGENOUS VARIABLES AND SHOCKS WITH SAME TIMING
% Recall that states are [exogenous(t);endogenous(t-1)]
PHI11=PHI2(1:nsh,1:nsh);
PHI21=PHI2(2*nsh+1:end,1:nsh);
PHI212=PHI2(2*nsh+1:end,nsh+1:2*nsh); %new gianni
PHI22=PHI2(2*nsh+1:end,2*nsh+1:end);

G1=[eye(nsh) zeros(nsh,neq-nsh);
    -PHI21 eye(neq-nsh)];
G2=[PHI11  zeros(nsh,neq-nsh);
   PHI212 PHI22];
Pi2(1:nsh,1:nsh)=eye(nsh);
RET=G1\G2;



%% Check the RE solution
res2=max(max(abs(A1*RET^2+A2*RET+A3))); % think of y(t-1)=Pi2*shocks; y(t)=Pi2*PHI2*shocks; y(t+1)=Pi2*PHI2^2*shocks



% LEARNING FROM THIS SOLUTION
% NOTE THAT Pi(1:nsh,:) is made of [I(:,end:-1:1),zeros] and the rest is
% identical to PHI(nsh+1:end,:). So need only learn about PHI2
% tmp=sum(abs(RET)>1e-10,1);
truenonzero0=find(sum(abs(RET)>=1e-10,1)>0);
truenonzero1=find(sum(abs(RET(truenonzero0,truenonzero0))>=1e-10,1)>0);
truenonzero0=truenonzero0(truenonzero1);
% truenonzero1=find(sum(abs(RET(:,truenonzero0))>1e-10,1)>0);
%              find(sum(abs(RET(:,truenonzero0(truenonzero1)))>1e-10,1)>0)
if ~isempty(const_st)
truenonzero=setdiff(truenonzero0,const_st);
else
    truenonzero=truenonzero0;
end
noshocks=input('Include shocks in the PLM?[1/0(default)] ==> ');
if isempty(noshocks)|| noshocks~=1
%%% IF SHOCKS ARE NOT OBSERVED ancomment the following
truenonzero=setdiff(truenonzero,1:nsh);
end
truezeros=setdiff(1:neq,truenonzero);

picknonzero=zeros(neq,1);
picknonzero(truenonzero)=1;
picknonzero=diag(picknonzero);
picknonzero=picknonzero(truenonzero,:);
picknojumps=zeros(neq,length(njump));
zz=1;
for kk=njump;
    picknojumps(kk,zz)=1;
    zz=zz+1;
end
% innovation loading matrix
Mshall=[Msh;zeros(neq-nsh,nsh)];

learning_sub;