Skip to content
Snippets Groups Projects
Commit 7ebc3985 authored by tzha's avatar tzha
Browse files

add two new m files

parent c2225930
Branches
No related tags found
No related merge requests found
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
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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment