1function [yhat,Estr,rcon,Ome,Rmat,u,v,d] = fidencond(valuecon,stepcon,varcon,nconstr,...
2                         nstepsm,nvar,lags,yfore,imf3s,phil,Bh,eq_iden,steps_iden)
3% Estimating conditional forecasting in the identified model
4%function [yhat,Estr,rcon,Ome,Rmat,u,v,d] = fidencond(valuecon,stepcon,varcon,nconstr,...
5%                         nstepsm,nvar,lags,yfore,imf3s,phil,Bh,eq_iden,steps_iden)
6%
7% valuecon:  vector of values conditioned
8% stepcon:   sequence (cell) of steps conditioned; if length(stepcon{i}) > 1, the condition
9%               is then an arithmetic average of log(y) over the stepcon{i} period.
10% varcon:    vector of variables conditioned
11% nconstr:   number of constraints
12% nstepsm:   maximum number of steps in all constraints
13% nvar:   number of variables in the BVAR model
14% lags:   number of lags in the BVAR model
15% yfore:  uncondtional forecasts: forep-by-nvar
16% imf3s: 3-dimensional impulse responses matrix:
17%                               impsteps-by-nvar shocks-by-nvar responses
18% phil:  the 1-by-(nvar*lags+1) data matrix where k=nvar*lags+1
19%                 (last period plus lags before the beginning of forecast)
20% Bh:  reduced-form parameter matrix: k-by-nvar, y(t) = X(t)*Bh+e(t)
21%                    where X(t) is k-by-nvar and y(t) is 1-by-nvar
22% eq_iden: identified equation or shock (in terms of number).
23%          If eq_iden=[], then'fidencond' is, similar to RATS, to compute
24%             forecasts with *all* shocks.
25%          If length(eq_iden)=1, compute forecasts with only "MS" shocks.
26%          If eq_iden = [a b c], a is # of constraints for all shocks (<=nconstr),
27%             b is location of MS, and c is max steps for all shocks.
28%             My recall is that this option has never been tested 9/18/98
29% steps_iden:  a vector (set) of steps for identified shocks.  Valid only
30%              if length(eq_iden)=1. Note, length(steps_iden) must nconstr.
31% ------
32% yhat:  conditional forecasts: forep-by-nvar
33% Estr:  backed-out structural shocks (from N(0,1))
34% rcon:  vector - the difference between valuecon and log(yfore) (unconditional forecasts)
35% Rcon:  k-by-q (q constranits and k=nvar*max(nsteps)) so that
36%                        Rcon'*e = rcon where e is k-by-1, where k=nvar*nstepm
37% Rmat:  nstepsm-by-nvar (shocks), for only one constraint at a time. See Zha's
38%                  Forecast (1), pp.5-6
39% Ome:  k-by-k: covariance matrix of vectorized structural shocks vec(Estr)
40%
41% [u,d,v]:  svd(Rcon,0)
42%
43%% See Zha's note "Forecast (1)" pp.5-7, RATS manual (some errors in RATS), etc.
44%
45%% Some notations:  y(t+1) = y(t)B1 + e(t+1)inv(A0). e(t+1) is 1-by-n.
46%%    Let r(t+1)=e(t+1)inv(A0) + e(t+2)C + .... where inv(A0) is impulse
47%%          response at t=1, C at t=2, etc. The row of inv(A0) or C is
48%%          all responses to one shock.
49%%    Let r be q-by-1 (such as r(1) = r(t+1)
50%%                 = y(t+1) (constrained) - y(t+1) (forecast)).
51%%    Use impulse responses to find out R (k-by-q) where k=nvar*nsteps
52%%        where nsteps the largest constrained step.  The key of the program
53%%        is to creat R using impulse responses
54%%    Optimal solution for shock e where R'*e=r and e is k-by-1 is
55%%                 e = R*inv(R'*R)*r.
56%
57% NOTE: the code needs to be improved, 10/19/98 (use lzpaper/fcstidcnd.m for the time being).
58%
59%
60% Copyright (C) 1997-2012 Tao Zha
61%
62% This free software: you can redistribute it and/or modify
63% it under the terms of the GNU General Public License as published by
64% the Free Software Foundation, either version 3 of the License, or
65% (at your option) any later version.
66%
67% It is distributed in the hope that it will be useful,
68% but WITHOUT ANY WARRANTY; without even the implied warranty of
69% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
70% GNU General Public License for more details.
71%
72% If you did not received a copy of the GNU General Public License
73% with this software, see <http://www.gnu.org/licenses/>.
74%
75
76
77%IdenShock = ~isempty(eq_iden);   % if not empty, the shock is identified
78forep=size(yfore,1);
79impsteps=size(imf3s,1);
80if (forep<nstepsm) | (impsteps<nstepsm)
81	disp('Increase # of forecast or impulse steps!!')
82   disp('Or decrease # of constraints (nconstr) or constrained steps (stepcon(i))!!')
83	error('Maximum of conditional steps > # of forecast or impulse steps!!')
84	%warning
85	%return
86end
87if max(steps_iden) < nconstr
88	disp('Increase # of identified steps or decrease # of constraints')
89	error('length(steps_iden) > nconstr')
90end
91
92kfs = nvar*nstepsm;   % k -- fs: free shocks
93%*** initializing
94Rcon = zeros(kfs,nconstr);   % R: k-by-q
95Econ = zeros(kfs,1);      % E: k-by-1
96rcon = zeros(nconstr,1);   % r: q-by-1
97%rcon=valuecon-diag(yfore(stepcon,varcon));  % another way is to use "loop" below.
98A0in = reshape(imf3s(1,:,:),nvar,nvar);  % nvar shocks-by-nvar responses
99
100
101for i=1:nconstr
102   rcon(i)=valuecon(i)-mean(yfore(stepcon{i},varcon(i)));
103	%rcon(i)=valuecon(i)-sum(yfore(stepcon{i},varcon(i)));
104   Rmat = zeros(nstepsm,nvar);
105           % Rmat: row--nstepsm, column--nvar shocks (here all shocks except
106	        %     the identified one are set to zero) for a particular
107           %     endogenous variable 'varcon(i)'.  See Zha Forcast (1), pp.6-7
108	r2mat = zeros(nstepsm,1);   % simply one identified equation
109           % Must be here inside the loop because it's matrix of one column of Rcon
110   for j=1:length(stepcon{i})
111      if (length(eq_iden)==1)
112			r2mat(1:stepcon{i}(j)) = r2mat(1:stepcon{i}(j)) + ...
113                                      imf3s(stepcon{i}(j):-1:1,eq_iden,varcon(i));
114      elseif (length(eq_iden)>1)
115         if (i<=eq_iden(1))
116            Rmat(1:stepcon{i}(j),:) = Rmat(1:stepcon{i}(j),:) + ...
117			                             imf3s(stepcon{i}(j):-1:1,:,varcon(i));
118         else
119            Rmat(1:eq_iden(3),:) = Rmat(1:eq_iden(3),:) + ...
120                     imf3s(stepcon{i}(j):-1:stepcon{i}(j)-eq_iden(3)+1,:,varcon(i));
121
122            Rmat(eq_iden(3)+1:stepcon{i}(j),eq_iden(2)) = ...
123                         Rmat(eq_iden(3)+1:stepcon{i}(j),eq_iden(2)) + ...
124                         imf3s(stepcon{i}(j)-eq_iden(3):-1:1,eq_iden(2),varcon(i));
125         end
126		else
127			Rmat(1:stepcon{i}(j),:) = Rmat(1:stepcon{i}(j),:) + ...
128			                             imf3s(stepcon{i}(j):-1:1,:,varcon(i));
129	                % Rmat: row--nstepsm, column--nvar shocks (here all shocks are
130						 %     *not* set to zero) for a particular endogenous
131                   %     variable 'varcon(i)'.  See Zha Forcast (1), pp.6-7
132		end
133   end
134	%
135	if (length(eq_iden)==1)
136      Rmat(steps_iden,eq_iden) = r2mat(steps_iden);   % for only one constraint at a time
137	end
138	Rmat=Rmat/length(stepcon{i});    % <<>>
139	Rmatt = Rmat';   % Now, nvar-by-nstepsm. I think here is where RATS has an error
140							  % i.e. "OVERR" is not transposed when overlaid to "CAPR"
141	Rcon(:,i)=Rmatt(:);      % Rcon: k-by-q where q=nconstr
142end
143
144if nconstr
145   [u d v]=svd(Rcon,0); %trial.   Rcon: k-by-q; u: k-by-q
146	% rtr = Rcon'*Rcon; %trial
147	% rtrinv = inv(Rcon'*Rcon); %trial
148	vd=v.*(ones(size(v,2),1)*diag(d)'); %trial
149	dinv = 1./diag(d);    % inv(diag(d))
150	vdinv=v.*(ones(size(v,2),1)*dinv'); %trial
151	rtr=vd*vd';       % R'*R
152	rtrinv = vdinv*vdinv';   % inv(R'*R)
153   Ome = eye(kfs) - u*u';       % note: I-u*u' = I - R*inv(R'*R)*R'; k-by-k
154
155	Econ = Rcon*rtrinv*rcon;    % E = R*inv(R'R)*r; mean
156else
157	Econ = zeros(kfs,1);
158	Rcon = NaN;
159	rcon = NaN;
160	u = NaN;
161	d = NaN;
162	v = NaN;
163end
164
165
166Estr = reshape(Econ,nvar,nstepsm);
167Estr = Estr';   % transpose so that
168          % Estr: structural shocks. Row--steps, Column--n shocks
169Ures = Estr*A0in;     % nstepsm-by-nvar
170			 % Ures: reduced-form residuals.  Row--steps; Column--n shocks
171%Ures = Estr*A0in;
172
173% ** reconstruct x(t) for y(t+h) = x(t+h-1)*B
174% **       where phi = x(t+h-1) with last column being constant
175%
176tcwc = nvar*lags;     % total coefficients without constant
177phi=phil;
178%
179yhat = zeros(forep,nvar);
180for k=1:forep
181	if (k<=nstepsm)
182   	epsl = Ures(k,:);
183   	yhat(k,:) = phi*Bh + epsl;
184		%yhat(k,:) = phi*Bh;
185   	phi(1,nvar+1:tcwc) = phi(1,1:tcwc-nvar);
186   	phi(1,1:nvar) = yhat(k,:);
187	else
188   	yhat(k,:) = phi*Bh;
189   	phi(1,nvar+1:tcwc) = phi(1,1:tcwc-nvar);
190   	phi(1,1:nvar) = yhat(k,:);
191	end
192end
193