diff --git a/MatlabFiles/MSV/fwz_msv_msre.m b/MatlabFiles/MSV/fwz_msv_msre.m
new file mode 100755
index 0000000000000000000000000000000000000000..79ffb6b25c74f535404e0b92889cab710a58c2ac
--- /dev/null
+++ b/MatlabFiles/MSV/fwz_msv_msre.m
@@ -0,0 +1,134 @@
+function [F1, F2, G1, G2, V, err] = fwz_msv_msre(P, A, B, Psi, s, x, max_count, tol)
+%[ F1 F2 G1 G2 V err ] = fwz_msv_msre(P, A, B, Psi, s, x, max_count, tol)
+% Computes MSV solution of
+%
+% A(s(t) x(t) = B(s(t) x(t-1) + Psi(s(t)) epsilon(t) + Pi eta(t)
+%
+% using Newton's Method. Assumes that Pi' = [zeros(s,n-s) eye(s)] and
+% that A{i} is invertible. P is the transition matrix and P(i,j) is the
+% probability that s(t+1)=j given that s(t)=i. Note that the rows of P
+% must sum to one. x is the initial value and if not passed is set to
+% zero. max_count is the maximum number of iterations on Newton's method
+% before failing and tol is the convergence criterion.
+%
+% The solution is of the form
+%
+% x(t) = V{s(t)}*F1{s(t)}*x(t-1) + V{s(t)}*G1{s(t)}*epsilon(t)
+%
+% eta(t) = F2{s(t)}*x(t-1) + G2{s(t)}*epsilont(t)
+%
+% err is set to 0 upon sucess. A positive value is the number of
+% iterations before the method was terminated.
+%
+
+h=size(P,1);
+n=size(A{1},1);
+r=size(Psi{1},2);
+
+if (nargin <= 7) || (tol <= 0)
+ tol=1e-5;
+end
+
+if (nargin <= 6) || (max_count <= 0)
+ max_count=1000;
+end
+
+if nargin <= 5
+ x=zeros(h*s*(n-s),1);
+end
+f=zeros(h*s*(n-s),1);
+
+
+Is=eye(s);
+I=eye(n-s);
+
+% The following code would be used if Pi were passed instead of the
+% assumption that Pi' = [zeros(s,n-s) eye(s)].
+% for i=1:h
+% [Q,R] = qr(Pi{i});
+% W=[zeros(n-s,s) I; inv(R(1:s,1:s)); zeros(s,n-s)]*Q';
+% A{i}=W*A{i};
+% B{i}=W*B{i};
+% Psi{i}=W*Psi{i};
+% end
+
+U=cell(h,1);
+for i=1:h
+ U{i}=inv(A{i});
+end
+C=cell(h,h);
+for i=1:h
+ for j=1:h
+ C{i,j}=P(i,j)*B{j}*U{i};
+ end
+end
+
+cont=true;
+count=1;
+D=zeros(h*s*(n-s),h*s*(n-s));
+X=cell(h,1);
+for i=1:h
+ X{i}=reshape(x((i-1)*s*(n-s)+1:i*s*(n-s)),s,n-s);
+end
+while cont
+ for i=1:h
+ for j=1:h
+ W1=C{i,j}*[I; -X{i}];
+ W2=W1(1:n-s,:);
+ D((i-1)*s*(n-s)+1:i*s*(n-s),(j-1)*s*(n-s)+1:j*s*(n-s)) = kron(W2',Is);
+ if i == j
+ W1=zeros(s,n);
+ for k=1:h
+ W1=W1+[X{k} Is]*C{i,k};
+ end
+ W2=-W1(:,n-s+1:end);
+ D((i-1)*s*(n-s)+1:i*s*(n-s),(j-1)*s*(n-s)+1:j*s*(n-s)) = D((i-1)*s*(n-s)+1:i*s*(n-s),(j-1)*s*(n-s)+1:j*s*(n-s)) + kron(I,W2);
+ end
+ end
+ end
+
+ for i=1:h
+ mf=zeros(s,n-s);
+ for j=1:h
+ mf=mf+[X{j} Is]*C{i,j}*[I; -X{i}];
+ end
+ f((i-1)*s*(n-s)+1:i*s*(n-s))=reshape(mf,s*(n-s),1);
+ end
+
+ y=D\f;
+ x=x - y;
+
+ if (count > max_count) || (norm(f) < tol)
+ cont=false;
+ else
+ count=count+1;
+ for i=1:h
+ X{i}=reshape(x((i-1)*s*(n-s)+1:i*s*(n-s)),s,n-s);
+ end
+ end
+end
+
+if (norm(f) < tol)
+ err=0;
+else
+ err=count;
+end
+
+F1=cell(h,1);
+F2=cell(h,1);
+G1=cell(h,1);
+G2=cell(h,1);
+V=cell(h,1);
+pi=[zeros(n-s,s); Is];
+for i=1:h
+ X=reshape(x((i-1)*s*(n-s)+1:i*s*(n-s)),s,n-s);
+ V{i}=U{i}*[I; -X];
+ W=[A{i}*V{i} pi];
+ F=W\B{i};
+ F1{i}=F(1:n-s,:);
+ F2{i}=F(n-s+1:end,:);
+ G=W\Psi{i};
+ G1{i}=G(1:n-s,:);
+ G2{i}=G(n-s+1:end,:);
+end
+
\ No newline at end of file
diff --git a/MatlabFiles/fn_multigraphn_3g_shadedbands.m b/MatlabFiles/fn_multigraphn_3g_shadedbands.m
new file mode 100755
index 0000000000000000000000000000000000000000..4a68388d53c4d78eaf15d38277e3b328a3a86d5d
--- /dev/null
+++ b/MatlabFiles/fn_multigraphn_3g_shadedbands.m
@@ -0,0 +1,171 @@
+function scaleout = fn_multigraphn_3g_shadedbands(imfn,...
+ xlab,ylab,tstring,XTick,YTickIndx,scaleIndx,nrowg,ncolg)
+%Searching <<>> for ad hoc and specific changes.
+%
+%Stacking n sets of impulse responses in one graph. See fn_multigraph2_ver2.m.
+% imfn: row: "nstp" time horizon (in the graphics),
+% column: "nrow "variables such as responses (row in the graphics),
+% 3rd D: across "ncol" different situations such as shocks (column in the graphics),
+% 4th D: low band, estimate, high band (in each graph).
+% xlab: x-axis labels on the top
+% ylab: y-axis labels on the left
+% tstring: string for time (e.g., month or quarter) -- x-axis labels at the bottom
+% YTickIndx: 1: enable YTick; 0: disable;
+% To get a better picture, it is sometimes better to set YtickIndx to zero.
+% scaleIndx: 1: enable scale along Y-axis; 0: disable
+% nrowg: number of rows in the graph
+% ncolg: number of columns in the graph
+%
+% See imrgraph, imcerrgraph, imrerrgraph
+nstp = size(imfn,1);
+nrow = size(imfn,2);
+ncol = size(imfn,3);
+nmodels = size(imfn,4);
+t = 1:nstp;
+treverse = fliplr(t);
+
+if (nargin < 9)
+ nrowg = nrow;
+ ncolg = ncol;
+end
+if (nrowg*ncolg < nrow*ncol)
+ nrowg
+ ncolg
+ nrow
+ ncol
+
+ error('fn_multigraphn_3g_shadedbands.m: nrowg*ncolg must be greater than nrow*ncol')
+end
+
+
+tempmax=zeros(ncol,1);
+tempmin=zeros(ncol,1);
+maxval=zeros(nrow,1);
+minval=zeros(nrow,1);
+for i = 1:nrow
+ for j = 1:ncol
+ tempmax(j) = -realmax;
+ tempmin(j) = realmax;
+ for k=1:nmodels
+ jnk = max(imfn(:,i,j,k));
+ tempmax(j) = max([jnk tempmax(j)]);
+ %
+ jnk = min(imfn(:,i,j,k));
+ tempmin(j) = min([jnk tempmin(j)]);
+ end
+ end
+ maxval(i)=max(tempmax);
+ minval(i)=min(tempmin);
+end
+
+scaleout = [maxval(:) minval(:)];
+
+%--------------
+% Column j: Shock 1 to N; Row i: Responses to
+%-------------
+%figure
+
+
+rowlabel = 1;
+for i = 1:nrow % column: from top to bottom
+ columnlabel = 1;
+
+ if minval(i)<0
+ if maxval(i)<=0
+ yt=[minval(i) 0];
+ else
+ yt=[minval(i) 0 maxval(i)];
+ end
+ else % (minval(i) >=0)
+ if maxval(i) > 0
+ yt=[0 maxval(i)];
+ else % (identically zero responses)
+ yt=[-1 0 1];
+ end
+ end
+
+
+ scale=[1 nstp minval(i) maxval(i)];
+
+ for j = 1:ncol % row: from left to right
+ k1=(i-1)*ncol+j;
+ subplot(nrowg,ncolg,k1)
+
+ %set(0,'DefaultAxesColorOrder',[1 0 0;0 1 0;0 0 1],...
+ % 'DefaultAxesLineStyleOrder','-|--|:')
+ %set(0,'DefaultAxesLineStyleOrder','-|--|:|-.')
+ %---<<>>
+ %set(0,'DefaultAxesColorOrder',[0 0 0],...
+ % 'DefaultAxesLineStyleOrder','-.|-.|-|--|-.|-*|--o|:d')
+
+ nseries = zeros(nstp, nmodels);
+ for k=1:nmodels
+ nseries(:,k) = imfn(:,i,j,k);
+ end
+
+ %---<<>>
+ fill([t treverse],[nseries(:,[3])' fliplr(nseries(:,[1])')],[0.8,0.8,0.8],'EdgeColor','none');
+ %plot(t,nseries(:,[1 3]), '-.k','LineWidth',1.0); %,'Color','k');
+ hold on
+ plot(t,nseries(:,[2]), '-k','LineWidth',1.5);
+ %plot(t,nseries(:,[4]), '--k','LineWidth',1.6);
+ hold off
+
+ %set(gca,'LineStyleOrder','-|--|:|-.')
+ %set(gca,'LineStyleOrder',{'-*',':','o'})
+ grid;
+ if scaleIndx
+ axis(scale);
+ end
+ %set(gca,'YLim',[minval(i) maxval(i)])
+ %
+ set(gca,'XTick',XTick)
+ if YTickIndx
+ set(gca,'YTick',yt)
+ end
+
+ if i<nrow
+ set(gca,'XTickLabelMode','manual','XTickLabel',[])
+ end
+ %set(gca,'XTickLabel',' ');
+ if (scaleIndx) && (j>1)
+ set(gca,'YTickLabel',' ');
+ end
+ if rowlabel == 1
+ %title(['x' num2str(j)])
+ %title(eval(['x' num2str(j)]))
+ if (~isempty(xlab)), title(char(xlab(j))), end
+ end
+ if columnlabel == 1
+ %ylabel(['x' num2str(i)])
+ %ylabel(eval(['x' num2str(i)]))
+ if (~isempty(ylab)), ylabel(char(ylab(i))), end
+ end
+ if (i==nrow) && (~isempty(tstring))
+ xlabel(tstring)
+ end
+ columnlabel = 0;
+ end
+ rowlabel = 0;
+end
+
+
+%Order of line styles and markers used in a plot.
+%This property specifies which line styles and markers to use and in what order
+%when creating multiple-line plots. For example,set(gca,'LineStyleOrder', '-*|:|o')sets LineStyleOrder to solid line with asterisk
+%marker, dotted line, and hollow circle marker. The default is (-), which specifies
+%a solid line for all data plotted. Alternatively, you can create a cell array
+%of character strings to define the line styles:set(gca,'LineStyleOrder',{'-*',':','o'})MATLAB supports four line styles, which you can specify any number of
+%times in any order. MATLAB cycles through the line styles only after using
+%all colors defined by the ColorOrder property. For example,
+%the first eight lines plotted use the different colors defined by ColorOrder with
+%the first line style. MATLAB then cycles through the colors again, using the
+%second line style specified, and so on.You can also specify line style and color directly with the plot and plot3 functions
+%or by altering the properties of theline or
+%lineseries objects after creating the graph. High-Level Functions and LineStyleOrderNote that, if the axes NextPlot property is set
+%to replace (the default), high-level functions like plot reset
+%the LineStyleOrder property before determining the line
+%style to use. If you want MATLAB to use a LineStyleOrder that
+%is different from the default, set NextPlot to replacechildren. Specifying a Default LineStyleOrderYou can also specify your own default LineStyleOrder.
+%For example, this statementset(0,'DefaultAxesLineStyleOrder',{'-*',':','o'})
+%creates a default value for