ms_mardd.m 8.23 KB
Newer Older
1
function ms_mardd(options_)
michel's avatar
michel committed
2
3
4
5
6
7
8
% Applies to both linear and exclusion restrictions.
% (1) Marginal likelihood function p(Y) for constant structural VAR models, using Chib (1995)'s ``Marginal Likelihood from the Gibbs Output'' in JASA.
% (2) Conditional likelihood function f(Y|A0, A+) on the ML estimate for constant exclusion-identified models.
% See Forecast (II) pp.67-80.
%
% Tao Zha, September 1999.  Quick revisions, May 2003.  Final revision, September 2004.

9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
% Copyright (C) 2011 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.

michel's avatar
michel committed
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
msstart2     % start the program in which everyhting is initialized through msstart2.m
if ~options_.ms.indxestima
   warning(' ')
   disp('You must set IxEstima=1 in msstart to run this program')
	disp('Press ctrl-c to abort now')
	pause
end

A0xhat = zeros(size(A0hat));
Apxhat = zeros(size(Aphat));
if (0)
   %Robustness check to see if the same result is obtained with the purterbation of the parameters.
   for k=1:nvar
      bk = Uiconst{k}'*A0hat(:,k);
      gk = Viconst{k}'*Aphat(:,k);
      A0xhat(:,k) =  Uiconst{k}*(bk + 5.2*randn(size(bk)));   % Perturbing the posterior estimate.
      Apxhat(:,k) = Viconst{k}*(gk + 5.2*randn(size(gk)));      % Perturbing the posterior estimate.
   end
else
   %At the posterior estimate.
   A0xhat = A0hat;   % ML estimate of A0
   Apxhat = Aphat;      % ML estimate of A+
end
%--- Rename variables.
YatYa = yty;
XatYa = xty;
ytx = xty';
YatXa = ytx;
XatXa = xtx;



%--------- The log value of p(A0,A+) at some point such as the peak ----------
vlog_a0p = 0;
Yexpt=0;  % exponential term for Y in p(Y|A0,A+) at some point such as the peak
Apexpt=0.0;  % 0.0 because we have chosen posterior estimate of A+ as A+*.  Exponential term for A+ conditional on A0 and Y
%======= Computing the log prior pdf of a0a+ and the exponential term for Y in p(Y|A0,A+).
for k=1:nvar
   a0k = A0xhat(:,k);  % meaningful parameters in the kth equation.
   apk = Apxhat(:,k);  % meaningful parameters in the kth equation.

   %--- Prior settings.
   S0bar = H0invtld{k};  %See Claim 2 on p.69b.
   Spbar = Hpinvtld{k};
   bk = Uiconst{k}'*a0k;  % free parameters in the kth equation.
   gk = Viconst{k}'*apk;  % free parameters in the kth equation.
   gbark = Ptld{k}*bk;   % bar: prior

   %--- The exponential term for Y in p(Y|A0,A+)
   Yexpt = Yexpt - 0.5*(a0k'*YatYa*a0k - 2*apk'*XatYa*a0k + apk'*XatXa*apk);
   %--- The log prior pdf.
   vlog_a0p = vlog_a0p - 0.5*(size(Uiconst{k},2)+size(Viconst{k},2))*log(2*pi) + 0.5*log(abs(det(S0bar))) + ...
         0.5*log(abs(det(Spbar))) - 0.5*(bk'*S0bar*bk+(gk-gbark)'*Spbar*(gk-gbark));
   %--- For p(A+|Y,a0) only.
   tmpd = gk - Pmat{k}*bk;
   Apexpt = Apexpt - 0.5*tmpd'*(Hpinv{k}*tmpd);
end
vlog_a0p

%--------- The log value of p(Y|A0,A+) at some point such as the peak. ----------
%--------- Note that logMarLHres is the same as vlog_Y_a, just to double check. ----------
vlog_Y_a = -0.5*nvar*fss*log(2*pi) + fss*log(abs(det(A0xhat))) + Yexpt
                 % a: given a0 and a+
logMarLHres = 0;   % Initialize log of the marginal likelihood (restricted or constant parameters).
for ki=1:fss   %ndobs+1:fss     % Forward recursion to get the marginal likelihood.  See F on p.19 and pp. 48-49.
   %----  Restricted log marginal likelihood function (constant parameters).
   [A0l,A0u] = lu(A0xhat);
   ada = sum(log(abs(diag(A0u))));   % log|A0|
   termexp = y(ki,:)*A0xhat - phi(ki,:)*Apxhat;   % 1-by-nvar
   logMarLHres = logMarLHres - (0.5*nvar)*log(2*pi) + ada - 0.5*termexp*termexp';  % log value
end
logMarLHres


%--------- The log value of p(A+|Y,A0) at some point such as the peak ----------
totparsp = 0.0;
tmpd = 0.0;
for k=1:nvar
   totparsp = totparsp + size(Viconst{k},2);
   tmpd = tmpd + 0.5*log(abs(det(Hpinv{k})));
end
vlog_ap_Ya0 = -0.5*totparsp*log(2*pi) + tmpd + Apexpt;




%===================================
%  Compute p(a0,k|Y,ao) at some point such as the peak (in this situation, we simply
%   generate results from the original Gibbs sampler).  See FORECAST (2) pp.70-71
%===================================
%--- Global set up for Gibbs.
[Tinv,UT] = fn_gibbsrvar_setup(H0inv, Uiconst, Hpinv, Pmat, Viconst, nvar, fss);
%
vlog_a0_Yao = zeros(nvar,1);
  % the log value of p(a0k|Y,ao) where ao: other a's at some point such as the peak of ONLY some a0's
vlog=zeros(ndraws2,1);
for k=1:nvar
   bk = Uiconst{k}'*A0xhat(:,k);
   indx_ks=[k:nvar];  % the columns that exclude 1-(k-1)th columns
   A0gbs0 = A0hat;   % starting at some point such as the peak
   nk = n0(k);

   if k<nvar
      %--------- The 1st set of draws to be tossed away. ------------------
      for draws = 1:ndraws1
         if ~mod(draws,nbuffer)
            disp(' ')
            disp(sprintf('The %dth column or equation in A0 with %d 1st tossed-away draws in Gibbs',k,draws))
         end
         A0gbs1 = fn_gibbsrvar(A0gbs0,UT,nvar,fss,n0,indx_ks);
         A0gbs0=A0gbs1;    % repeat the Gibbs sampling
      end


      %--------- The 2nd set of draws to be used. ------------------
      for draws = 1:ndraws2
         if ~mod(draws,nbuffer)
            disp(' ')
            disp(sprintf('The %dth column or equation in A0 with %d usable draws in Gibbs',k,draws))
         end
         [A0gbs1, Wcell] = fn_gibbsrvar(A0gbs0,UT,nvar,fss,n0,indx_ks);
         %------ See p.71, Forecast (II).
         %------ Computing p(a0_k|Y,a_others) at some point such as the peak along the dimensions of indx_ks.
         Vk = Tinv{k}\Wcell{k};  %V_k on p.71 of Forecast (II).
         gbeta = Vk\bk;  % inv(V_k)*b_k on p.71 of Forecast (II) where alpha_k = b_k in our notation.
         [Vtq,Vtr]=qr(Vk',0);  %To get inv(V_k)'*inv(V_k) in (*) on p.71 of Forecast (II).
         %
         vlog(draws) = 0.5*(fss+nk)*log(fss)-log(abs(det(Vk)))-0.5*(nk-1)*log(2*pi)-...
                  0.5*(fss+1)*log(2)-gammaln(0.5*(fss+1))+fss*log(abs(gbeta(1)))-...
                  0.5*fss*bk'*(Vtr\(Vtr'\bk));

         A0gbs0=A0gbs1;    % repeat the Gibbs sampling
      end
      vlogm=max(vlog);
      qlog=vlog-vlogm;
      vlogxhat=vlogm-log(ndraws2)+log(sum(exp(qlog)));
      vlog_a0_Yao(k) = vlogxhat;
         % The log value of p(a0_k|Y,a_others) where a_others: other a's at some point such as the peak of ONLY some a0's
   else
      disp(' ')
      disp(sprintf('The last(6th) column or equation in A0 with no Gibbs draws'))
      [A0gbs1, Wcell] = fn_gibbsrvar(A0gbs0,UT,nvar,fss,n0,indx_ks)
      %------ See p.71, Forecast (II).
      %------ Computing p(a0_k|Y,a_others) at some point such as the peak along the dimensions of indx_ks.
      Vk = Tinv{k}\Wcell{k};  %V_k on p.71 of Forecast (II).
      gbeta = Vk\bk;  % inv(V_k)*b_k on p.71 of Forecast (II) where alpha_k = b_k in our notation.
      [Vtq,Vtr]=qr(Vk',0);  %To get inv(V_k)'*inv(V_k) in (*) on p.71 of Forecast (II).
      %
      vloglast = 0.5*(fss+nk)*log(fss)-log(abs(det(Vk)))-0.5*(nk-1)*log(2*pi)-...
               0.5*(fss+1)*log(2)-gammaln(0.5*(fss+1))+fss*log(abs(gbeta(1)))-...
               0.5*fss*bk'*(Vtr\(Vtr'\bk));
      vlog_a0_Yao(k) = vloglast;
   end
end
ndraws2

disp('Prior pdf -- log(p(a0hat, a+hat)):');
vlog_a0p
disp('LH pdf -- log(p(Y|a0hat, a+hat)):');
vlog_Y_a
disp('Posterior Kernal -- logp(ahat) + logp(Y|ahat):');
vlog_Y_a + vlog_a0p
disp('Posterior pdf -- log(p(a0_i_hat|a0_other_hat, Y)):');
vlog_a0_Yao
disp('Posterior pdf -- log(p(aphat|a0hat, Y)):');
vlog_ap_Ya0

%--------- The value of marginal density p(Y) ----------
disp(' ');
disp(' ');
disp('************ Marginal Likelihood of Y or Marginal Data Density: ************');
vlogY = vlog_a0p+vlog_Y_a-sum(vlog_a0_Yao)-vlog_ap_Ya0