diff --git a/.gitignore b/.gitignore
index b102fd00ea9bf051d74d31206f78f4934508dbbe..eba503751ecd72c5482cd06f6b81aedd10949b8c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -21,6 +21,7 @@
*.npy
*.o
*.dll
+*.zip
# Can also ignore all directories and files in a directory.
# tmp/**/*
diff --git a/MatlabFiles/MSV/OldVersions/dw_gensys.m b/MatlabFiles/MSV/OldVersions/dw_gensys.m
new file mode 100755
index 0000000000000000000000000000000000000000..cd959321a68ebe04c981adf3bdaac052548533f2
--- /dev/null
+++ b/MatlabFiles/MSV/OldVersions/dw_gensys.m
@@ -0,0 +1,46 @@
+function [G1,impact,gev,eu]=dw_gensys(g0,g1,psi,pi,div)
+% function [G1,impact,gev,eu]=gensys(g0,g1,c,psi,pi,div)
+% System given as
+% g0*y(t)=g1*y(t-1)+psi*e(t)+pi*eta(t),
+% with e an exogenous variable process and eta being endogenously determined
+% one-step-ahead expectational errors. Returned system is
+% y(t)=G1*y(t-1)+impact*e(t).
+% If div is omitted from argument list, a div>1 is calculated.
+% eu(1)=1 for existence, eu(2)=1 for uniqueness. eu(1)=-1 for
+% existence only with not-s.c. z; eu=[-2,-2] for coincident zeros.
+% By Daniel Waggoner
+
+realsmall=1e-8;
+
+G1=[];
+impact=[];
+gev=[];
+
+n=size(pi,1);
+
+if nargin == 4
+ div=1.0;
+end
+
+[a b q z]=qz(g0,g1);
+
+for i=1:n
+ if (abs(a(i,i)) < realsmall) & (abs(b(i,i)) < realsmall)
+ disp('Coincident zeros.')
+ eu=[-2;-2];
+ gev=[diag(a) diag(b)];
+ return
+ end
+end
+
+[a b q z]=qzsort(a,b,q,z);
+
+suppress=0;
+i=n;
+while (i > 0) & (abs(b(i,i)) > abs(a(i,i)*div))
+ suppress=suppress+1;
+ i=i-1;
+end
+
+[G1,impact,eu]=qzcomputesolution(a,b,q,z,psi,pi,suppress);
+
diff --git a/MatlabFiles/MSV/OldVersions/gensys.m b/MatlabFiles/MSV/OldVersions/gensys.m
new file mode 100755
index 0000000000000000000000000000000000000000..69f7fc0fa0229de91d3d4cc9a2be063d24152c42
--- /dev/null
+++ b/MatlabFiles/MSV/OldVersions/gensys.m
@@ -0,0 +1,214 @@
+function [G1,C,impact,fmat,fwt,ywt,gev,eu]=gensys(g0,g1,c,psi,pi,div)
+% function [G1,C,impact,fmat,fwt,ywt,gev,eu]=gensys(g0,g1,c,psi,pi,div)
+% System given as
+% g0*y(t)=g1*y(t-1)+c+psi*z(t)+pi*eta(t),
+% with z an exogenous variable process and eta being endogenously determined
+% one-step-ahead expectational errors. Returned system is
+% y(t)=G1*y(t-1)+C+impact*z(t)+ywt*inv(I-fmat*inv(L))*fwt*z(t+1) .
+% If z(t) is i.i.d., the last term drops out.
+% If div is omitted from argument list, a div>1 is calculated.
+% eu(1)=1 for existence, eu(2)=1 for uniqueness. eu(1)=-1 for
+% existence only with not-s.c. z; eu=[-2,-2] for coincident zeros.
+% By Christopher A. Sims
+% Corrected 10/28/96 by CAS
+
+eu=[0;0];
+realsmall=1e-7;
+fixdiv=(nargin==6);
+n=size(g0,1);
+[a b q z v]=qz(g0,g1);
+if ~fixdiv, div=1.01; end
+nunstab=0;
+zxz=0;
+for i=1:n
+% ------------------div calc------------
+ if ~fixdiv
+ if abs(a(i,i)) > 0
+ divhat=abs(b(i,i))/abs(a(i,i));
+ % bug detected by Vasco Curdia and Daria Finocchiaro, 2/25/2004 A root of
+ % exactly 1.01 and no root between 1 and 1.02, led to div being stuck at 1.01
+ % and the 1.01 root being misclassified as stable. Changing < to <= below fixes this.
+ if 1+realsmall<divhat & divhat<=div
+ div=.5*(1+divhat);
+ end
+ end
+ end
+% ----------------------------------------
+ nunstab=nunstab+(abs(b(i,i))>div*abs(a(i,i)));
+ if abs(a(i,i))<realsmall & abs(b(i,i))<realsmall
+ zxz=1;
+ end
+end
+div ;
+nunstab;
+if ~zxz
+ [a b q z]=qzdiv(div,a,b,q,z);
+end
+gev=[diag(a) diag(b)];
+if zxz
+ disp('Coincident zeros. Indeterminacy and/or nonexistence.')
+ eu=[-2;-2];
+ % correction added 7/29/2003. Otherwise the failure to set output
+ % arguments leads to an error message and no output (including eu).
+ G1=[];C=[];impact=[];fmat=[];fwt=[];ywt=[];gev=[];
+ return
+end
+q1=q(1:n-nunstab,:);
+q2=q(n-nunstab+1:n,:);
+z1=z(:,1:n-nunstab)';
+z2=z(:,n-nunstab+1:n)';
+a2=a(n-nunstab+1:n,n-nunstab+1:n);
+b2=b(n-nunstab+1:n,n-nunstab+1:n);
+etawt=q2*pi;
+% zwt=q2*psi;
+[ueta,deta,veta]=svd(etawt);
+md=min(size(deta));
+bigev=find(diag(deta(1:md,1:md))>realsmall);
+ueta=ueta(:,bigev);
+veta=veta(:,bigev);
+deta=deta(bigev,bigev);
+% ------ corrected code, 3/10/04
+eu(1) = length(bigev)>=nunstab;
+length(bigev);
+nunstab;
+% ------ Code below allowed "existence" in cases where the initial lagged state was free to take on values
+% ------ inconsistent with existence, so long as the state could w.p.1 remain consistent with a stable solution
+% ------ if its initial lagged value was consistent with a stable solution. This is a mistake, though perhaps there
+% ------ are situations where we would like to know that this "existence for restricted initial state" situation holds.
+%% [uz,dz,vz]=svd(zwt);
+%% md=min(size(dz));
+%% bigev=find(diag(dz(1:md,1:md))>realsmall);
+%% uz=uz(:,bigev);
+%% vz=vz(:,bigev);
+%% dz=dz(bigev,bigev);
+%% if isempty(bigev)
+%% exist=1;
+%% else
+%% exist=norm(uz-ueta*ueta'*uz) < realsmall*n;
+%% end
+%% if ~isempty(bigev)
+%% zwtx0=b2\zwt;
+%% zwtx=zwtx0;
+%% M=b2\a2;
+%% for i=2:nunstab
+%% zwtx=[M*zwtx zwtx0];
+%% end
+%% zwtx=b2*zwtx;
+%% [ux,dx,vx]=svd(zwtx);
+%% md=min(size(dx));
+%% bigev=find(diag(dx(1:md,1:md))>realsmall);
+%% ux=ux(:,bigev);
+%% vx=vx(:,bigev);
+%% dx=dx(bigev,bigev);
+%% existx=norm(ux-ueta*ueta'*ux) < realsmall*n;
+%% else
+%% existx=1;
+%% end
+% ----------------------------------------------------
+% Note that existence and uniqueness are not just matters of comparing
+% numbers of roots and numbers of endogenous errors. These counts are
+% reported below because usually they point to the source of the problem.
+% ------------------------------------------------------
+[ueta1,deta1,veta1]=svd(q1*pi);
+md=min(size(deta1));
+bigev=find(diag(deta1(1:md,1:md))>realsmall);
+ueta1=ueta1(:,bigev);
+veta1=veta1(:,bigev);
+deta1=deta1(bigev,bigev);
+%% if existx | nunstab==0
+%% %disp('solution exists');
+%% eu(1)=1;
+%% else
+%% if exist
+%% %disp('solution exists for unforecastable z only');
+%% eu(1)=-1;
+%% %else
+%% %fprintf(1,'No solution. %d unstable roots. %d endog errors.\n',nunstab,size(ueta1,2));
+%% end
+%% %disp('Generalized eigenvalues')
+%% %disp(gev);
+%% %md=abs(diag(a))>realsmall;
+%% %ev=diag(md.*diag(a)+(1-md).*diag(b))\ev;
+%% %disp(ev)
+%% % return;
+%% end
+if isempty(veta1)
+ unique=1;
+else
+ unique=norm(veta1-veta*veta'*veta1)<realsmall*n;
+end
+
+veta1;
+veta;
+norm(veta1-veta*veta'*veta1);
+
+if unique
+ %disp('solution unique');
+ eu(2)=1;
+else
+ fprintf(1,'Indeterminacy. %d loose endog errors.\n',size(veta1,2)-size(veta,2));
+ %disp('Generalized eigenvalues')
+ %disp(gev);
+ %md=abs(diag(a))>realsmall;
+ %ev=diag(md.*diag(a)+(1-md).*diag(b))\ev;
+ %disp(ev)
+% return;
+end
+tmat = [eye(n-nunstab) -(ueta*(deta\veta')*veta1*deta1*ueta1')'];
+G0= [tmat*a; zeros(nunstab,n-nunstab) eye(nunstab)];
+G1= [tmat*b; zeros(nunstab,n)];
+% ----------------------
+% G0 is always non-singular because by construction there are no zeros on
+% the diagonal of a(1:n-nunstab,1:n-nunstab), which forms G0's ul corner.
+% -----------------------
+G0I=inv(G0);
+G1=G0I*G1;
+usix=n-nunstab+1:n;
+C=G0I*[tmat*q*c;(a(usix,usix)-b(usix,usix))\q2*c];
+impact=G0I*[tmat*q*psi;zeros(nunstab,size(psi,2))];
+fmat=b(usix,usix)\a(usix,usix);
+fwt=-b(usix,usix)\q2*psi;
+ywt=G0I(:,usix);
+% -------------------- above are output for system in terms of z'y -------
+G1=real(z*G1*z');
+C=real(z*C);
+impact=real(z*impact);
+% Correction 10/28/96: formerly line below had real(z*ywt) on rhs, an error.
+ywt=z*ywt;
+
+%--------------------------------------------------------------------------
+% The code below check that a solution is obtained with stable manifold
+% z1. Should be commented out for production runs.
+%--------------------------------------------------------------------------
+% eu
+%
+% z2=z(:,n-nunstab+1:n);
+% z1=null(z2');
+%
+% n=norm(g0*(G1*z1) - g1*z1);
+% disp('norm(g0*(g1*z1) - g1*z1) -- should be zero');
+% disp(n);
+%
+% eta=-veta*diag(1./diag(deta))*ueta'*q2*psi;
+% n=norm(g0*impact - (psi + pi*eta));
+% disp('norm(g0*impact - psi - pi*eta) -- should be zero');
+% disp(n);
+%
+% n=norm(z2'*G1*z1);
+% disp('norm(z2''*G1*z1) -- should be zero');
+% disp(n);
+%
+% n=norm(z2'*G1);
+% disp('norm(z2''G1) -- best if zero');
+% disp(n);
+%
+% n=norm(z2'*impact);
+% disp('norm(z2''*impact) -- should be zero');
+% disp(n);
+%
+% n=norm(G1*z2);
+% disp('norm(G1*z2)');
+% disp(n);
+%
+% disp('press any key to continue');
+% pause
diff --git a/MatlabFiles/MSV/OldVersions/mainfunction.m b/MatlabFiles/MSV/OldVersions/mainfunction.m
new file mode 100755
index 0000000000000000000000000000000000000000..cf38e917c721f806fccd86e33fe94d6bb181ee4e
--- /dev/null
+++ b/MatlabFiles/MSV/OldVersions/mainfunction.m
@@ -0,0 +1,219 @@
+% The constant paramater model written in the canonical form:
+% A x_t = B x_{t-1} + Gamma_u u_t + Pi pi_t
+% where f_t = [pi_t, y_t, R_t, mu_w_t, nu_t, E_t pi_{t+1}, E_t y_{t+1}]';
+% ut = [Sst*x_t{-1}; diag(e_t(1),...,e_t(h*))]
+% e_t is a 3-by-1 vector of fundamental shocks: e_{r t}, e_{wt}, and e_{nu t};
+% pi_t is a 2-by-1 vector of expectationsl errors for inflation and output: pi_{pi t} and pi_{y t}.
+%
+% Output format: x_t = G1*x_{t-1} + impact*e_t.
+%
+% Look for <<>> for ad hoc changes.
+%See Notes pp.41-42.
+
+smallval = 1.0e-8; %<sqrt(machine epsilon).
+locs_normalizedvars = [3 4 5]; %R_t, mu_w_t, nu_t;
+seednumber = 7910; %472534; % if 0, random state at each clock time
+if seednumber
+ randn('state',seednumber);
+ rand('state',seednumber);
+else
+ randn('state',fix(100*sum(clock)));
+ rand('state',fix(100*sum(clock)));
+end
+n_normalizedvars = length(locs_normalizedvars);
+
+
+
+%--- Reading in the values of the parameters.
+%addpath D:\ZhaData\WorkDisk\LiuWZ\CalibrationMatlab\LWZmodel\Paravals2r
+cd ./Paravals2r
+%paravals2r_bench
+%paravals2r_complex
+paravals2r_noexistence
+%paravals2r_try
+cd ..
+
+Q = [
+ q11 1-q22
+ 1-q11 q22
+ ];
+disp(' ')
+disp('Ergodic probability:')
+qbar = fn_ergodp(Q)
+g_etabar = qbar'*g_eta;
+if (flag_swap_rs)
+ g_eta = flipud(g_eta);
+ g_gamma = flipud(g_gamma);
+ g_rhor = flipud(g_rhor);
+ g_phipi = flipud(g_phipi);
+ g_phiy = flipud(g_phiy);
+ Q = flipud(fliplr(Q));
+end
+
+
+%--- Composite regime-switching parameters. The symbol "td" stands for "tilde."
+htd = 1; %4;
+Qtd = [
+ Q(1,1) Q(1,1) 0 0
+ 0 0 Q(1,2) Q(1,2)
+ Q(2,1) Q(2,1) 0 0
+ 0 0 Q(2,2) Q(2,2)
+ ];
+%--- psi1(std_t) = g_psi1(s_t,s_{t-1})
+psi1td = [
+ g_etabar*(1-g_eta(1)) / (g_eta(1)*(1-g_eta(1)))
+ g_etabar*(1-g_eta(2)) / (g_eta(2)*(1-g_eta(1)))
+ g_etabar*(1-g_eta(1)) / (g_eta(1)*(1-g_eta(2)))
+ g_etabar*(1-g_eta(2)) / (g_eta(2)*(1-g_eta(2)))
+ ];
+%--- psi2(std_t) = g_psi2(s_{t-1})
+psi2td = [
+ (1-g_beta*g_etabar)*(1-g_eta(1)) / (g_eta(1)*(1+g_thetap*(1-g_alpha)/g_alpha));
+ (1-g_beta*g_etabar)*(1-g_eta(2)) / (g_eta(2)*(1+g_thetap*(1-g_alpha)/g_alpha));
+ (1-g_beta*g_etabar)*(1-g_eta(1)) / (g_eta(1)*(1+g_thetap*(1-g_alpha)/g_alpha));
+ (1-g_beta*g_etabar)*(1-g_eta(2)) / (g_eta(2)*(1+g_thetap*(1-g_alpha)/g_alpha));
+ ];
+%--- gamma0(std_t) = g_gamma(s_{t})
+gamma0td = [
+ g_gamma(1)
+ g_gamma(1)
+ g_gamma(2)
+ g_gamma(2)
+ ];
+%--- gamma1(std_t) = g_gamma(s_{t-1})
+gamma1td = [
+ g_gamma(1)
+ g_gamma(2)
+ g_gamma(1)
+ g_gamma(2)
+ ];
+%--- rhor(std_t) = g_rhor(s_{t})
+rhortd = [
+ g_rhor(1)
+ g_rhor(1)
+ g_rhor(2)
+ g_rhor(2)
+ ];
+%--- phipi(std_t) = g_phipi(s_{t})
+phipitd = [
+ g_phipi(1)
+ g_phipi(1)
+ g_phipi(2)
+ g_phipi(2)
+ ];
+%--- phiy(std_t) = g_phiy(s_{t})
+phiytd = [
+ g_phiy(1)
+ g_phiy(1)
+ g_phiy(2)
+ g_phiy(2)
+ ];
+
+%-------------------- Filling in A(std_t)'s, B(std_t)'s, etc for an expanded gensys. -------------------%
+Astd = cell(htd,1);
+Bstd = cell(htd,1);
+Psistd = cell(htd,1);
+for si=1:htd
+ a1td = (1+g_beta*psi1td(si)*gamma0td(si));
+ a2td = psi2td(si)*((g_xi+1)/g_alpha + b/(g_lambda-b));
+ Astd{si} = [
+ -a1td a2td 0 psi2td(si) psi2td(si)*(b/(g_lambda-b)) g_beta*psi1td(si) 0
+ 0 -(g_lambda+b)/g_lambda -(g_lambda-b)/g_lambda 0 g_rhov - b/g_lambda (g_lambda-b)/g_lambda 1
+ -(1-rhortd(si))*phipitd(si) -(1-rhortd(si))*phiytd(si) 1 0 0 0 0
+ 0 0 0 1 0 0 0
+ 0 0 0 0 1 0 0
+ ];
+ Bstd{si} = [
+ -gamma1td(si) psi2td(si)*b/(g_lambda-b) 0 0 0 0 0
+ 0 -b/g_lambda 0 0 0 0 0
+ 0 0 rhortd(si) 0 0 0 0
+ 0 0 0 g_rhow 0 0 0
+ 0 0 0 0 g_rhov 0 0
+ ];
+ Psistd{si} = [
+ 0 0 0
+ 0 0 0
+ g_sigmar 0 0
+ 0 g_sigmaw 0
+ 0 0 g_sigmav
+ ];
+end
+%--- <<>>Reordering
+if (flag_swap_fps)
+ ki = 4;
+ kk = 1;
+ Atmp = Astd{ki};
+ Astd{ki} = Astd{kk};
+ Astd{kk} = Atmp;
+ Btmp = Bstd{ki};
+ Bstd{ki} = Bstd{kk};
+ Bstd{kk} = Btmp;
+ SwapM = zeros(ki,ki);
+ SwapM(kk,ki) = kk;
+ SwapM(ki,1) = 1;
+ SwapM(2,2) = 1;
+ SwapM(3,3) = 1;
+ Qtd = SwapM*Qtd*SwapM;
+end
+
+
+
+
+Aind = cell(htd,1);
+Bind = cell(htd,1);
+Psiind = cell(htd,1);
+Piind = cell(htd,1);
+G1ind = cell(htd,1);
+impactind = cell(htd,1);
+gevind = cell(htd,1);
+gevind_abs = cell(htd,1);
+euind = cell(htd,1);
+for si=1:htd
+ Aind{si} = [
+ Astd{si}
+ 1 0 0 0 0 0 0
+ 0 1 0 0 0 0 0
+ ];
+ Bind{si} = [
+ Bstd{si}
+ 0 0 0 0 0 1 0
+ 0 0 0 0 0 0 1
+ ];
+ Psiind{si} = [
+ Psistd{si}
+ 0 0 0
+ 0 0 0
+ ];
+ Piind{si} = [
+ zeros(5,2)
+ 1 0
+ 0 1
+ ];
+
+ disp('---------- Individual (possibly swapped) regimes in isolation --------------')
+ disp('Under Regime')
+ si
+
+ %[G1ind{si},Cind,impactind{si},fmat,fwt,ywt,gevind{si},euind{si}]=gensys(Aind{si}, Bind{si}, zeros(7,1), Psiind{si}, Piind{si}, divind(si));
+ %[G1ind{si},impactind{si},gevind{si},euind{si}]=msv_one(Aind{si}, Bind{si}, Psiind{si}, Piind{si});
+ [G1ind{si},impactind{si},gevind{si},euind{si}]=msv_simple(Aind{si}, Bind{si}, Psiind{si}, Piind{si});
+ [G1 impact gev eu]=msv_all(Aind{si}, Bind{si}, Psiind{si}, Piind{si});
+ %[G1ind{si},impactind{si},v2junk,gevind{si},euind{si}]=gensys_dw(Aind{si}, Bind{si}, Psiind{si}, Piind{si}, divind(si));
+
+ k=size(eu,2);
+ for i=1:k
+ gev_abs = abs(gev{i}(:,1).\gev{i}(:,2));
+ gev_abs
+ G1{i}
+ impact{i}
+ eu{i}
+ end
+end
+% disp(' ')
+% disp('----------------- Existence and uniqueness for each individual regime -----------------------')
+% for si=1:htd
+% euind{si}
+% end
+return
+
+
diff --git a/MatlabFiles/MSV/OldVersions/msv_all.m b/MatlabFiles/MSV/OldVersions/msv_all.m
new file mode 100755
index 0000000000000000000000000000000000000000..7fffafc11a21b995b35c8d85ba81484aca939a72
--- /dev/null
+++ b/MatlabFiles/MSV/OldVersions/msv_all.m
@@ -0,0 +1,142 @@
+function [G1,impact,gev,eu]=msv_all(g0,g1,psi,pi)
+% function [G1,C,impact,gev,eu]=msv_all(g0,g1,psi,pi)
+% System given as
+% g0*y(t)=g1*y(t-1)+psi*z(t)+pi*eta(t),
+% with z an exogenous variable process and eta being endogenously determined
+% one-step-ahead expectational errors. Returned system is
+% y(t)=G1*y(t-1)+impact*z(t).
+% eu(1)=1 for existence, eu(2)=1 for uniqueness.
+% eu=[-2,-2] for coincident zeros.
+% Attempts all possible ordering such that all explosive roots are
+% suppressed and no complex conjugate pairs are split
+% By Daniel Waggoner -- Based on code by Christopher A. Sims
+
+realsmall=1e-10;
+
+n=size(pi,1);
+s=size(pi,2);
+ii=1;
+
+[a b q z]=qz(g0,g1);
+
+for i=1:n
+ if (abs(a(i,i)) < realsmall) & (abs(b(i,i)) < realsmall)
+ disp('Coincident zeros.')
+ eu{ii}=[-2;-2];
+ gev{ii}=[diag(a) diag(b)];
+ G1{ii}=[]
+ impact{ii}=[]
+ return
+ end
+end
+
+[a b q z]=qzsort(a,b,q,z);
+
+% determine complex conjugate pairs
+pairs=zeros(n,1);
+i=1;
+while i < n
+ if abs(b(i+1,i+1)*conj(a(i,i)) - a(i+1,i+1)*conj(b(i,i))) < realsmall
+ pairs(i)=1;
+ pairs(i+1)=-1;
+ i=i+2;
+ else
+ i=i+1;
+ end
+end
+
+% determine number unstable roots
+nunstable=0;
+i=1;
+while i <= n
+ if pairs(i) == 0
+ if (abs(b(i,i)) > abs(a(i,i)))
+ nunstable=n-i+1;
+ break;
+ else
+ i=i+1;
+ end
+ else
+ if (abs(b(i,i)) > abs(a(i,i))) | (abs(b(i+1,i+1)) > abs(a(i+1,i+1)))
+ nunstable=n-i+1;
+ break;
+ else
+ i=i+2;
+ end
+ end
+end
+
+% No bounded MSV solutions?
+if nunstable > s
+ disp('No bounded MSV solutions')
+ eu{ii}=[0;0];
+ gev{ii}=[diag(a) diag(b)];
+ G1{ii}=[]
+ impact{ii}=[]
+ return;
+end
+
+% Multiple bounded MSV solutions?
+if nunstable < s
+ disp('Multiple bounded MSV solutions should exist');
+
+ % Get roots to suppress
+ idx=zeros(n,1);
+ idx_size=zeros(n,1);
+ k=1;
+ idx(1)=n-nunstable;
+ cont=1;
+ while cont > 0
+ if (pairs(idx(k)) == -1)
+ idx_size(k)=2;
+ else
+ idx_size(k)=1;
+ end
+
+ d=nunstable;
+ for i=1:k
+ d=d+idx_size(i);
+ end
+
+ if (d == s)
+ cont=cont+1;
+ [an bn qn zn]=qzmoveindex(a,b,q,z,pairs,idx,k,n-nunstable);
+ [G1{ii},impact{ii},eu{ii}]=qzcomputesolution(an,bn,qn,zn,psi,pi,s);
+ gev{ii}=[diag(an) diag(bn)];
+ ii=ii+1;
+ end
+
+ if (d >= s)
+ if idx(k) > idx_size(k)
+ idx(k)=idx(k)-idx_size(k);
+ else
+ if k > 1
+ k=k-1;
+ idx(k)=idx(k)-idx_size(k);
+ else
+ return;
+ end
+ end
+ else
+ if idx(k) > idx_size(k)
+ idx(k+1)=idx(k)-idx_size(k);
+ k=k+1;
+ else
+ if k > 1
+ k=k-1;
+ idx(k)=idx(k)-idx_size(k);
+ else
+ return;
+ end
+ end
+ end
+ end
+end
+
+% Unique MSV solution?
+disp('Unique bounded MVS solution');
+[G1{ii},impact{ii},eu{ii}]=qzcomputesolution(a,b,q,z,psi,pi,s);
+gev{ii}=[diag(a) diag(b)];
+
+
+
diff --git a/MatlabFiles/MSV/OldVersions/msv_one.m b/MatlabFiles/MSV/OldVersions/msv_one.m
new file mode 100755
index 0000000000000000000000000000000000000000..44c79456919075887ebfc5a122a53f6ff84d6711
--- /dev/null
+++ b/MatlabFiles/MSV/OldVersions/msv_one.m
@@ -0,0 +1,145 @@
+function [G1,impact,gev,eu]=msv_one(g0,g1,psi,pi)
+% function [G1,C,impact,gev,eu]=msv_one(g0,g1,psi,pi)
+% System given as
+% g0*y(t)=g1*y(t-1)+psi*z(t)+pi*eta(t),
+% with z an exogenous variable process and eta being endogenously determined
+% one-step-ahead expectational errors. Returned system is
+% y(t)=G1*y(t-1)+impact*z(t).
+% eu(1)=1 for existence, eu(2)=1 for uniqueness. eu(1)=-1 for
+% existence only with not-s.c. z; eu=[-2,-2] for coincident zeros.
+% By Daniel Waggoner -- Based on code by Christopher A. Sims
+
+realsmall=1e-10;
+
+n=size(pi,1);
+s=size(pi,2);
+
+[a b q z]=qz(g0,g1);
+
+for i=1:n
+ if (abs(a(i,i)) < realsmall) & (abs(b(i,i)) < realsmall)
+ disp('Coincident zeros.')
+ eu=[-2;-2];
+ gev=[diag(a) diag(b)];
+ G1=[]
+ impact=[]
+ return
+ end
+end
+
+[a b q z]=qzsort(a,b,q,z);
+
+% determine complex conjugate pairs
+pairs=zeros(n,1);
+i=1;
+while i < n
+ if abs(b(i+1,i+1)*conj(a(i,i)) - a(i+1,i+1)*conj(b(i,i))) < realsmall
+ pairs(i)=1;
+ pairs(i+1)=-1;
+ i=i+2;
+ else
+ i=i+1;
+ end
+end
+
+% determine number unstable roots
+nunstable=0;
+i=1;
+while i <= n
+ if pairs(i) == 0
+ if (abs(b(i,i)) > abs(a(i,i)))
+ nunstable=n-i+1;
+ break;
+ else
+ i=i+1;
+ end
+ else
+ if (abs(b(i,i)) > abs(a(i,i))) | (abs(b(i+1,i+1)) > abs(a(i+1,i+1)))
+ nunstable=n-i+1;
+ break;
+ else
+ i=i+2;
+ end
+ end
+end
+
+% No bounded MSV solutions?
+if nunstable > s
+ disp('No bounded MSV solutions')
+ eu=[0;0];
+ gev=[diag(a) diag(b)];
+ G1=[]
+ impact=[]
+ return;
+end
+
+% Multiple bounded MSV solutions?
+if nunstable < s
+ disp('Multiple bounded MSV solutions should exist');
+
+ % Get roots to suppress
+ idx=zeros(n,1);
+ idx_size=zeros(n,1);
+ k=1;
+ idx(1)=n-nunstable;
+ cont=1;
+ while cont > 0
+ if (pairs(idx(k)) == -1)
+ idx_size(k)=2;
+ else
+ idx_size(k)=1;
+ end
+
+ d=nunstable;
+ for i=1:k
+ d=d+idx_size(i);
+ end
+
+ if (d == s)
+ cont=cont+1;
+ [an bn qn zn]=qzmoveindex(a,b,q,z,pairs,idx,k,n-nunstable);
+ [G1,impact,eu]=qzcomputesolution(an,bn,qn,zn,psi,pi,s);
+ if (eu(1) == 1) & (eu(2) == 1)
+ gev=[diag(an) diag(bn)];
+ return;
+ else
+ disp('MSV solution does not exist for the ordering:');
+ abs(diag(bn)./diag(an))
+ eu
+ end
+ end
+
+ if (d >= s)
+ if idx(k) > idx_size(k)
+ idx(k)=idx(k)-idx_size(k);
+ else
+ if k > 1
+ k=k-1;
+ idx(k)=idx(k)-idx_size(k);
+ else
+ disp('No MSV solutions exist');
+ return;
+ end
+ end
+ else
+ if idx(k) > idx_size(k)
+ idx(k+1)=idx(k)-idx_size(k);
+ k=k+1;
+ else
+ if k > 1
+ k=k-1;
+ idx(k)=idx(k)-idx_size(k);
+ else
+ disp('No MSV solutions exist');
+ return;
+ end
+ end
+ end
+ end
+end
+
+% Unique MSV solution?
+disp('Unique bounded MVS solution');
+[G1,impact,eu]=qzcomputesolution(a,b,q,z,psi,pi,s);
+gev=[diag(a) diag(b)];
+
diff --git a/MatlabFiles/MSV/OldVersions/msv_simple.m b/MatlabFiles/MSV/OldVersions/msv_simple.m
new file mode 100755
index 0000000000000000000000000000000000000000..c37148fc3e6e279cf95ee4e33a669458e45cc6a4
--- /dev/null
+++ b/MatlabFiles/MSV/OldVersions/msv_simple.m
@@ -0,0 +1,62 @@
+function [G1,impact,gev,eu]=msv_simple(g0,g1,psi,pi)
+% function [G1,C,impact,gev,eu]=msv_one(g0,g1,psi,pi)
+% System given as
+% g0*y(t)=g1*y(t-1)+psi*z(t)+pi*eta(t),
+% with z an exogenous variable process and eta being endogenously determined
+% one-step-ahead expectational errors. Returned system is
+% y(t)=G1*y(t-1)+impact*z(t).
+% eu(1)=1 for existence, eu(2)=1 for uniqueness.
+% eu=[-2,-2] for coincident zeros.
+% By Daniel Waggoner -- Based on code by Christopher A. Sims
+
+realsmall=1e-10;
+
+G1=[];
+impact=[];
+gev=[];
+
+n=size(pi,1);
+s=size(pi,2);
+
+[a b q z]=qz(g0,g1);
+
+for i=1:n
+ if (abs(a(i,i)) < realsmall) & (abs(b(i,i)) < realsmall)
+ disp('Coincident zeros.')
+ eu=[-2;-2];
+ gev=[diag(a) diag(b)];
+ return
+ end
+end
+
+[a b q z]=qzsort(a,b,q,z);
+
+% determine number unstable roots
+nunstable=0;
+i=n;
+while i > 0
+ if (abs(b(i,i)) > abs(a(i,i)))
+ nunstable=nunstable+1;
+ i=i-1;
+ else
+ break;
+ end
+end
+
+eu=[0;0];
+gev=[diag(a) diag(b)];
+
+% No bounded MSV solutions?
+if nunstable > s
+ disp('No bounded MSV solutions')
+ return;
+end
+
+% Multiple or unique bounded MSV solutions?
+if nunstable < s
+ disp('Multiple bounded MSV solutions should exist');
+else
+ disp('Unique bounded MVS solution');
+end
+
+[G1,impact,eu]=qzcomputesolution(a,b,q,z,psi,pi,s);
diff --git a/MatlabFiles/MSV/OldVersions/qzcomputesolution.m b/MatlabFiles/MSV/OldVersions/qzcomputesolution.m
new file mode 100755
index 0000000000000000000000000000000000000000..2312735e29e3b833ffd66d2129635022e0f5eaf9
--- /dev/null
+++ b/MatlabFiles/MSV/OldVersions/qzcomputesolution.m
@@ -0,0 +1,126 @@
+function [G1,impact,eu]=qzcomputesolution(a,b,q,z,psi,pi,suppress)
+% function [G1,C,impact,eu]=qzcomputesolution(a,b,q,z,c,psi,pi,suppress)
+% System given as
+% (q'*a*z')*y(t)=(q'*b*z')*y(t-1)+psi*e(t)+pi*eta(t),
+% with e an exogenous variable process and eta an endogenously determined
+% one-step-ahead expectational error. q and z are unitary matrices and a
+% and b are upper triangular matrices that are computed via a generalized
+% Schur decomposition. It is assumed that the ordering of the system
+% is such that complex conjugate pairs are together. If a complex
+% conjugate pair is split, then G1 and impact will be complex.
+% Returned system is
+% y(t)=G1*y(t-1)+impact*z(t).
+% eu(1)=1 for existence, eu(2)=1 for uniqueness.
+%
+% The last suppress roots are suppressed.
+%
+% By Daniel Waggoner -- Based on code by Christopher A. Sims
+
+realsmall=sqrt(eps);
+
+eu=[0;0];
+G1=[];
+impact=[];
+
+n=size(pi,1);
+s=size(pi,2);
+
+q1=q(1:n-suppress,:);
+q2=q(n-suppress+1:n,:);
+
+[u1 d1 v1]=svd(q2*pi,'econ');
+if suppress > 0
+ m=sum(diag(d1) > realsmall*d1(1,1));
+else
+ m=0;
+end
+u1=u1(:,1:m);
+d1=d1(1:m,1:m);
+v1=v1(:,1:m);
+
+if (m == s)
+ eu(2)=1;
+else
+ eu(2)=0;
+end
+
+[u2 d2 v2]=svd(q2*psi,'econ');
+if suppress > 0
+ m=sum(diag(d2) > realsmall*d2(1,1));
+else
+ m=0;
+end
+u2=u2(:,1:m);
+d2=d2(1:m,1:m);
+v2=v2(:,1:m);
+
+if norm((eye(suppress) - u1*u1')*u2) < realsmall
+ eu(1)=1;
+else
+ eu(1)=0;
+ return;
+end
+
+
+% Compute impact and G0
+a11_inv=inv(a(1:n-suppress,1:n-suppress));
+x=a11_inv*q1*pi*v1*diag(1./diag(d1))*u1';
+impact=[a11_inv*q1*psi - x*q2*psi; zeros(suppress,size(psi,2))];
+
+G1=zeros(n,n);
+G1(1:n-suppress,1:n-suppress)=a11_inv*b(1:n-suppress,1:n-suppress);
+
+% this is to make the answer agree with gensys - this piece does not effect
+% the answer, at least after the initial period.
+G1(1:n-suppress,n-suppress+1:n)=a11_inv*b(1:n-suppress,n-suppress+1:n) - x*b(n-suppress+1:n,n-suppress+1:n);
+
+% Convert back to y
+impact=z*impact;
+G1=z*G1*z';
+
+% Is the solution real?
+if (suppress == 0) | (suppress == n) | (abs(b(n-suppress+1,n-suppress+1)*conj(a(n-suppress,n-suppress)) - a(n-suppress+1,n-suppress+1)*conj(b(n-suppress,n-suppress))) > realsmall)
+ impact=real(impact);
+ G1=real(G1);
+end
+
+%--------------------------------------------------------------------------
+% The code below check that a solution is obtained with stable manifold
+% z1. Should be commented out for production runs.
+%--------------------------------------------------------------------------
+% z2=z(:,n-suppress+1:n);
+% z1=null(z2');
+% A=q'*a*z';
+% B=q'*b*z';
+%
+% n=norm(A*(G1*z1) - B*z1);
+% disp('norm(A*(g1*z1) - B*z1) -- should be zero');
+% disp(n);
+%
+% eta=-v1*diag(1./diag(d1))*u1'*q2*psi;
+% n=norm(A*impact - psi - pi*eta);
+% disp('norm(A*impact - psi - pi*eta) -- should be zero');
+% disp(n);
+%
+% n=norm(z2'*G1*z1);
+% disp('norm(z2''*G1*z1) -- should be zero');
+% disp(n);
+%
+% n=norm(z2'*G1);
+% disp('norm(z2''G1) -- best if zero');
+% disp(n);
+%
+% n=norm(z2'*impact);
+% disp('norm(z2''*impact) -- should be zero');
+% disp(n);
+%
+% n=norm(G1*z2);
+% disp('norm(G1*z2)');
+% disp(n);
+%
+% disp('press any key to continue');
+% pause
+
+
+
+
diff --git a/MatlabFiles/MSV/OldVersions/qzdiv.m b/MatlabFiles/MSV/OldVersions/qzdiv.m
new file mode 100755
index 0000000000000000000000000000000000000000..b04068fe95afabf3cf4a096e90849529fef410d9
--- /dev/null
+++ b/MatlabFiles/MSV/OldVersions/qzdiv.m
@@ -0,0 +1,40 @@
+function [A,B,Q,Z,v] = qzdiv(stake,A,B,Q,Z,v)
+%function [A,B,Q,Z,v] = qzdiv(stake,A,B,Q,Z,v)
+%
+% Takes U.T. matrices A, B, orthonormal matrices Q,Z, rearranges them
+% so that all cases of abs(B(i,i)/A(i,i))>stake are in lower right
+% corner, while preserving U.T. and orthonormal properties and Q'AZ' and
+% Q'BZ'. The columns of v are sorted correspondingly.
+%
+% by Christopher A. Sims
+% modified (to add v to input and output) 7/27/00
+vin = nargin==6;
+%if ~vin, v=[], end;
+if ~vin, v=[]; end;
+[n jnk] = size(A);
+root = abs([diag(A) diag(B)]);
+root(:,1) = root(:,1)-(root(:,1)<1.e-13).*(root(:,1)+root(:,2));
+root(:,2) = root(:,2)./root(:,1);
+for i = n:-1:1
+ m=0;
+ for j=i:-1:1
+ if (root(j,2) > stake | root(j,2) < -.1)
+ m=j;
+ break
+ end
+ end
+ if (m==0)
+ return
+ end
+ for k=m:1:i-1
+ [A B Q Z] = qzswitch(k,A,B,Q,Z);
+ tmp = root(k,2);
+ root(k,2) = root(k+1,2);
+ root(k+1,2) = tmp;
+ if vin
+ tmp=v(:,k);
+ v(:,k)=v(:,k+1);
+ v(:,k+1)=tmp;
+ end
+ end
+end
diff --git a/MatlabFiles/MSV/OldVersions/qzmoveindex.m b/MatlabFiles/MSV/OldVersions/qzmoveindex.m
new file mode 100755
index 0000000000000000000000000000000000000000..853b79d90263ac3d76dff67f466edab1f4f217a3
--- /dev/null
+++ b/MatlabFiles/MSV/OldVersions/qzmoveindex.m
@@ -0,0 +1,17 @@
+function [a,b,q,z] = qzmoveindex(a,b,q,z,pairs,idx,k,j)
+% function [a,b,q,z] = qzmoveindex(a,b,q,z,pairs,idx,k,j)
+%
+% Takes U.T. matrices a, b, orthonormal matrices q,z, rearranges them
+% so that the indices in idx are moved up to the jth position, while
+% preserving U.T. and orthogonal properties and q'az' and q'bz'.
+%
+% by Daniel Waggoner based on code by Christopher A. Sims
+
+for i=1:k
+ [a b q z pairs]=qzslide(a,b,q,z,pairs,idx(i),j);
+ if (pairs(j) == -1)
+ j=j-2;
+ else
+ j=j-1;
+ end
+end
\ No newline at end of file
diff --git a/MatlabFiles/MSV/OldVersions/qzslide.m b/MatlabFiles/MSV/OldVersions/qzslide.m
new file mode 100755
index 0000000000000000000000000000000000000000..f9bbcb46d93fd903c5f3e3aa31c81d9ce49c2f51
--- /dev/null
+++ b/MatlabFiles/MSV/OldVersions/qzslide.m
@@ -0,0 +1,80 @@
+function [a,b,q,z,pairs] = qzslide(a,b,q,z,pairs,i,j)
+% function [a,b,q,z,pairs] = qzslide(a,b,q,z,pairs,i,j)
+%
+% Takes U.T. matrices A, B, orthonormal matrices Q,Z, rearranges them
+% so that the ith diagonal element is moved to the jth position, while
+% preserving U.T. and orthogonal properties and Q'AZ' and Q'BZ'.
+%
+% The array pairs is also rearranged. Complex conjugate pairs are kept
+% adjacent. If pairs(i)=1, then position i and i+1 are complex
+% conjugate pairs. If pairs(i)=-1, then position i and i-1 are complex
+% conjugate pairs. If pairs(i)=0, then position i is real.
+%
+%
+% by Daniel Waggoner based on code by Christopher A. Sims
+
+n=size(a,1);
+
+if i == j
+ return;
+end
+
+if i < j
+ if pairs(j) == 1
+ j=j+1;
+ end
+ if pairs(i) == -1
+ j=j-1;
+ i=i-1;
+ else
+ if (pairs(i) == 1) & (j == n)
+ if (i == n-1)
+ return;
+ end
+ j=n-1;
+ end
+ end
+ while i < j
+ if (pairs(i) == 1)
+ [a b q z]=qzswitch(i+1,a,b,q,z);
+ [a b q z]=qzswitch(i,a,b,q,z);
+ pairs(i)=pairs(i+2);
+ pairs(i+1)=1;
+ pairs(i+2)=-1;
+ else
+ [a b q z]=qzswitch(i,a,b,q,z);
+ pairs(i)=pairs(i+1);
+ pairs(i+1)=0;
+ end
+ i=i+1;
+ end
+else
+ if pairs(j) == -1
+ j=j-1;
+ end
+ if pairs(i) == 1
+ j=j+1;
+ i=i+1;
+ else
+ if (pairs(i) == -1) & (j == 1)
+ if (i == 2)
+ return;
+ end
+ j=2;
+ end
+ end
+ while i > j
+ if (pairs(i) == -1)
+ [a b q z]=qzswitch(i-2,a,b,q,z);
+ [a b q z]=qzswitch(i-1,a,b,q,z);
+ pairs(i)=pairs(i-2);
+ pairs(i-2)=1;
+ pairs(i-1)=-1;
+ else
+ [a b q z]=qzswitch(i-1,a,b,q,z);
+ pairs(i)=pairs(i-1);
+ pairs(i-1)=0;
+ end
+ i=i-1;
+ end
+end
diff --git a/MatlabFiles/MSV/OldVersions/qzsort.m b/MatlabFiles/MSV/OldVersions/qzsort.m
new file mode 100755
index 0000000000000000000000000000000000000000..7f24331a9650d14c0ae42c25fb2f1a67bed4c754
--- /dev/null
+++ b/MatlabFiles/MSV/OldVersions/qzsort.m
@@ -0,0 +1,21 @@
+function [A,B,Q,Z] = qzsort(A,B,Q,Z)
+%function [A,B,Q,Z] = qzsort(stake,A,B,Q,Z)
+%
+% Takes U.T. matrices A, B, orthonormal matrices Q,Z, rearranges them
+% so that abs(B(i,i)/A(i,i)) are increasing, while preserving U.T. and
+% orthonormal properties and Q'AZ' and % Q'BZ'.
+%
+% by Daniel Waggoner based on code by Christopher A. Sims
+
+n=size(A,1);
+
+for i=2:n
+ for j=i:-1:2
+ if abs(A(j,j)*B(j-1,j-1)) > abs(A(j-1,j-1)*B(j,j))
+ [A B Q Z]=qzswitch(j-1,A,B,Q,Z);
+ else
+ break;
+ end
+ end
+end
+
diff --git a/MatlabFiles/MSV/OldVersions/qzswitch.m b/MatlabFiles/MSV/OldVersions/qzswitch.m
new file mode 100755
index 0000000000000000000000000000000000000000..93c36640cc3488a66dc6c5261b8778e45c314c94
--- /dev/null
+++ b/MatlabFiles/MSV/OldVersions/qzswitch.m
@@ -0,0 +1,60 @@
+function [A,B,Q,Z] = qzswitch(i,A,B,Q,Z)
+%function [A,B,Q,Z] = qzswitch(i,A,B,Q,Z)
+%
+% Takes U.T. matrices A, B, orthonormal matrices Q,Z, interchanges
+% diagonal elements i and i+1 of both A and B, while maintaining
+% Q'AZ' and Q'BZ' unchanged. If diagonal elements of A and B
+% are zero at matching positions, the returned A will have zeros at both
+% positions on the diagonal. This is natural behavior if this routine is used
+% to drive all zeros on the diagonal of A to the lower right, but in this case
+% the qz transformation is not unique and it is not possible simply to switch
+% the positions of the diagonal elements of both A and B.
+ realsmall=sqrt(eps)*10;
+%realsmall=1e-3;
+a = A(i,i); d = B(i,i); b = A(i,i+1); e = B(i,i+1);
+c = A(i+1,i+1); f = B(i+1,i+1);
+ % A(i:i+1,i:i+1)=[a b; 0 c];
+ % B(i:i+1,i:i+1)=[d e; 0 f];
+if (abs(c)<realsmall & abs(f)<realsmall)
+ if abs(a)<realsmall
+ % l.r. coincident 0's with u.l. of A=0; do nothing
+ return
+ else
+ % l.r. coincident zeros; put 0 in u.l. of a
+ wz=[b; -a];
+ wz=wz/sqrt(wz'*wz);
+ wz=[wz [wz(2)';-wz(1)'] ];
+ xy=eye(2);
+ end
+elseif (abs(a)<realsmall & abs(d)<realsmall)
+ if abs(c)<realsmall
+ % u.l. coincident zeros with l.r. of A=0; do nothing
+ return
+ else
+ % u.l. coincident zeros; put 0 in l.r. of A
+ wz=eye(2);
+ xy=[c -b];
+ xy=xy/sqrt(xy*xy');
+ xy=[[xy(2)' -xy(1)'];xy];
+ end
+else
+ % usual case
+ wz = [c*e-f*b, (c*d-f*a)'];
+ xy = [(b*d-e*a)', (c*d-f*a)'];
+ n = sqrt(wz*wz');
+ m = sqrt(xy*xy');
+ if m<eps*100
+ % all elements of A and B proportional
+ return
+ end
+ wz = n\wz;
+ xy = m\xy;
+ wz = [wz; -wz(2)', wz(1)'];
+ xy = [xy;-xy(2)', xy(1)'];
+end
+A(i:i+1,:) = xy*A(i:i+1,:);
+B(i:i+1,:) = xy*B(i:i+1,:);
+A(:,i:i+1) = A(:,i:i+1)*wz;
+B(:,i:i+1) = B(:,i:i+1)*wz;
+Z(:,i:i+1) = Z(:,i:i+1)*wz;
+Q(i:i+1,:) = xy*Q(i:i+1,:);
\ No newline at end of file
diff --git a/MatlabFiles/MSV/OldVersions/test_gensys.m b/MatlabFiles/MSV/OldVersions/test_gensys.m
new file mode 100755
index 0000000000000000000000000000000000000000..f8ae311cd2304f84cb4bedfc81871f8a161a837d
--- /dev/null
+++ b/MatlabFiles/MSV/OldVersions/test_gensys.m
@@ -0,0 +1,53 @@
+r=2;
+s=2;
+n=6;
+explosive=s;
+
+div=1.0;
+
+cont=1;
+while cont > 0
+ [q tmp]=qr(rand(n,n));
+ [z tmp]=qr(rand(n,n));
+ a=q*diag(ones(n,1)+rand(n,1))*z;
+ %a=q*z;
+ d=rand(n,1);
+ d(1:explosive)=d(1:explosive)+ones(explosive,1);
+ b=q*diag(d)*z;
+
+ psi=rand(n,r);
+ pi=rand(n,s);
+ c=zeros(n,1);
+
+ [G1,C,impact,fmat,fwt,ywt,gev,eu]=gensys(a,b,c,psi,pi,div);
+ [dwG1,dwimpact,dwgev,dweu]=dw_gensys(a,b,psi,pi,div);
+
+ G1
+ dwG1
+
+ impact
+ dwimpact
+
+ eu
+ dweu
+
+ if eu(1) == 1
+ cont=0;
+ end
+end
+
+% Simulate
+y0=rand(n,1)
+n_sim=10;
+for i=1:10
+ epsilon=rand(r,1);
+ y=G1*y0+impact*epsilon;
+ dwy=dwG1*y0+dwimpact*epsilon;
+
+ y
+ dwy
+ norm(y-dwy)
+ pause
+
+ y0=dwy;
+end
\ No newline at end of file
diff --git a/MatlabFiles/MSV/OldVersions/test_msv.m b/MatlabFiles/MSV/OldVersions/test_msv.m
new file mode 100755
index 0000000000000000000000000000000000000000..c408620adca36353d8b6b98e5b78fcf96b1fc4b5
--- /dev/null
+++ b/MatlabFiles/MSV/OldVersions/test_msv.m
@@ -0,0 +1,35 @@
+r=2;
+s=2;
+n=6;
+explosive=s;
+
+div=1.0;
+
+cont=1;
+while cont > 0
+ [q tmp]=qr(rand(n,n));
+ [z tmp]=qr(rand(n,n));
+ a=q*diag(ones(n,1)+rand(n,1))*z;
+ %a=q*z;
+ d=rand(n,1);
+ d(1:explosive)=d(1:explosive)+ones(explosive,1);
+ b=q*diag(d)*z;
+
+ psi=rand(n,r);
+ pi=rand(n,s);
+
+ [G1,G2,gev,eu]=msv_all(a,b,psi,pi);
+
+ k=size(eu,2);
+
+ for i=1:k
+ eu{i}
+ abs(gev{i}(:,2)./gev{i}(:,1))
+ G1{i}
+ G2{i}
+ i
+ k
+ pause
+ end
+end
+
diff --git a/MatlabFiles/MSV/msv_all_complex_AR.m b/MatlabFiles/MSV/msv_all_complex_AR.m
new file mode 100755
index 0000000000000000000000000000000000000000..3ad0abddcc6caed9624bbea5d2b41fce42827d6d
--- /dev/null
+++ b/MatlabFiles/MSV/msv_all_complex_AR.m
@@ -0,0 +1,112 @@
+function [G1,Impact,gev,z2,err]=msv_all_complex_AR(A,B,Psi,Pi)
+% System given as
+%
+% A*y(t) = B*y(t-1) + psi*epsilon(t) + pi*eta(t),
+%
+% with epsilon(t) an exogenous process and eta(t) endogenously determined
+% one-step-ahead expectational errors. Returned solutions are
+%
+% y(t)=G1{i}*y(t-1)+Impact{i}*epsilon(t).
+%
+% The span of the each of the solutions returned will be n-s and the
+% solution will be uniquely determined by its span. gev returns the
+% generalized eigenvalues, of which the last s were suppressed. z2 is the
+% subspace that is perpendicular to the span of the solution.
+%
+% err is the error return.
+% err = 0 -- success, at least one msv solution found.
+% err = 1 -- no msv solutions
+% err = 2 -- degenerate system
+%
+% Solutions are obtained via calls to qzcomputesolution().
+%
+% By Daniel Waggoner
+
+realsmall=sqrt(eps);
+
+n=size(Pi,1);
+s=size(Pi,2);
+
+[a,b,q,z]=qz(A,B,'complex');
+
+G1=cell(0,1);
+Impact=cell(1,1);
+gev=cell(0,1);
+gev{1}=[diag(a) diag(b)];
+z2=cell(0,1);
+err=1;
+
+% determine if the system is degenerate and count the number of diagonal
+% elements of a that are non-zero.
+kk=0;
+for i=1:n
+ if abs(a(i,i)) < realsmall
+ if abs(b(i,i)) < realsmall
+ disp('Coincident zeros.')
+ err=2;
+ return
+ end
+ else
+ kk=kk+1;
+ end
+end
+
+% too many roots must be suppressed. no msv type solutions.
+if (n-kk > s)
+ return
+end
+
+% sort the generalized eigenvalues
+[d,id]=sort(abs(diag(b)./diag(a)),'descend');
+[id,idx]=sort(id);
+[a0,b0,q0,z0]=ordqz(a,b,q,z,idx);
+a=a0; b=b0; q=q0; z=z0;
+
+idx=ones(n,1);
+for i=0:s-1
+ idx(n-i)=0;
+end
+
+ii=0;
+cont=1;
+cntmax = 500;
+cntnumber = 1;
+while (cont == 1) && (cntnumber <=cntmax)
+ % compute solution
+ [g1,impact,eu]=qzcomputesolution(a,b,q,z,Psi,Pi,s);
+
+ % save solution if unique
+ if (eu == [1;1])
+ ii=ii+1;
+ G1{ii,1}=g1;
+ Impact{ii,1}=impact;
+ gev{ii,1}=[diag(a) diag(b)];
+ z2{ii,1}=z(:,n-s+1:n);
+ err=0;
+ end
+
+ % increment idx
+ cont=0;
+ jj=1;
+ while (jj < kk) & (idx(jj) == 0)
+ idx(jj)=1;
+ jj=jj+1;
+ end
+ for j=jj+1:kk
+ if idx(j) == 0
+ idx(j)=1;
+ cont=1;
+ break;
+ end
+ end
+ for i=1:jj
+ idx(j-i)=0;
+ end
+
+ % reorder roots
+ [a,b,q,z]=ordqz(a0,b0,q0,z0,idx);
+
+ %--- Ad hoc addtion.
+ cntnumber = cntnumber+1;
+end
+
diff --git a/tzDocumentation/readme_use_dwswitch.prn b/tzDocumentation/readme_use_dwswitch.prn
index c62c1ccb1a06b31fc4c989862098cbf53ebce08a..42558ee0e36d4e6bf532a4331ce10b437daa8166 100755
--- a/tzDocumentation/readme_use_dwswitch.prn
+++ b/tzDocumentation/readme_use_dwswitch.prn
@@ -24,7 +24,7 @@
/* How to call base regimes
/* How to write out transition matrices
/*
-/* DSGE output files
+/* DSGE output files such as simulations and impulse responses
*****************************
@@ -684,7 +684,7 @@ fgetc(stdin); //fgetc(stdin) will automatically fflush(stdout). Use ctrl-c to
//==================================
-// DSGE output files
+// DSGE output files such as simulations and impulse responses
//==================================
Output format for ir_percentiles_TAG.prn reporting impulse responses with error bands
-------------------------------------------------------------------------------------------------------