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.
0001 function ms_mardd(options_) 0002 % Applies to both linear and exclusion restrictions. 0003 % (1) Marginal likelihood function p(Y) for constant structural VAR models, using Chib (1995)'s ``Marginal Likelihood from the Gibbs Output'' in JASA. 0004 % (2) Conditional likelihood function f(Y|A0, A+) on the ML estimate for constant exclusion-identified models. 0005 % See Forecast (II) pp.67-80. 0006 % 0007 % Tao Zha, September 1999. Quick revisions, May 2003. Final revision, September 2004. 0008 0009 % Copyright (C) 2011 Dynare Team 0010 % 0011 % This file is part of Dynare. 0012 % 0013 % Dynare is free software: you can redistribute it and/or modify 0014 % it under the terms of the GNU General Public License as published by 0015 % the Free Software Foundation, either version 3 of the License, or 0016 % (at your option) any later version. 0017 % 0018 % Dynare is distributed in the hope that it will be useful, 0019 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0020 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0021 % GNU General Public License for more details. 0022 % 0023 % You should have received a copy of the GNU General Public License 0024 % along with Dynare. If not, see <http://www.gnu.org/licenses/>. 0025 0026 msstart2 % start the program in which everyhting is initialized through msstart2.m 0027 if ~options_.ms.indxestima 0028 warning(' ') 0029 disp('You must set IxEstima=1 in msstart to run this program') 0030 disp('Press ctrl-c to abort now') 0031 pause 0032 end 0033 0034 A0xhat = zeros(size(A0hat)); 0035 Apxhat = zeros(size(Aphat)); 0036 if (0) 0037 %Robustness check to see if the same result is obtained with the purterbation of the parameters. 0038 for k=1:nvar 0039 bk = Uiconst{k}'*A0hat(:,k); 0040 gk = Viconst{k}'*Aphat(:,k); 0041 A0xhat(:,k) = Uiconst{k}*(bk + 5.2*randn(size(bk))); % Perturbing the posterior estimate. 0042 Apxhat(:,k) = Viconst{k}*(gk + 5.2*randn(size(gk))); % Perturbing the posterior estimate. 0043 end 0044 else 0045 %At the posterior estimate. 0046 A0xhat = A0hat; % ML estimate of A0 0047 Apxhat = Aphat; % ML estimate of A+ 0048 end 0049 %--- Rename variables. 0050 YatYa = yty; 0051 XatYa = xty; 0052 ytx = xty'; 0053 YatXa = ytx; 0054 XatXa = xtx; 0055 0056 0057 0058 %--------- The log value of p(A0,A+) at some point such as the peak ---------- 0059 vlog_a0p = 0; 0060 Yexpt=0; % exponential term for Y in p(Y|A0,A+) at some point such as the peak 0061 Apexpt=0.0; % 0.0 because we have chosen posterior estimate of A+ as A+*. Exponential term for A+ conditional on A0 and Y 0062 %======= Computing the log prior pdf of a0a+ and the exponential term for Y in p(Y|A0,A+). 0063 for k=1:nvar 0064 a0k = A0xhat(:,k); % meaningful parameters in the kth equation. 0065 apk = Apxhat(:,k); % meaningful parameters in the kth equation. 0066 0067 %--- Prior settings. 0068 S0bar = H0invtld{k}; %See Claim 2 on p.69b. 0069 Spbar = Hpinvtld{k}; 0070 bk = Uiconst{k}'*a0k; % free parameters in the kth equation. 0071 gk = Viconst{k}'*apk; % free parameters in the kth equation. 0072 gbark = Ptld{k}*bk; % bar: prior 0073 0074 %--- The exponential term for Y in p(Y|A0,A+) 0075 Yexpt = Yexpt - 0.5*(a0k'*YatYa*a0k - 2*apk'*XatYa*a0k + apk'*XatXa*apk); 0076 %--- The log prior pdf. 0077 vlog_a0p = vlog_a0p - 0.5*(size(Uiconst{k},2)+size(Viconst{k},2))*log(2*pi) + 0.5*log(abs(det(S0bar))) + ... 0078 0.5*log(abs(det(Spbar))) - 0.5*(bk'*S0bar*bk+(gk-gbark)'*Spbar*(gk-gbark)); 0079 %--- For p(A+|Y,a0) only. 0080 tmpd = gk - Pmat{k}*bk; 0081 Apexpt = Apexpt - 0.5*tmpd'*(Hpinv{k}*tmpd); 0082 end 0083 vlog_a0p 0084 0085 %--------- The log value of p(Y|A0,A+) at some point such as the peak. ---------- 0086 %--------- Note that logMarLHres is the same as vlog_Y_a, just to double check. ---------- 0087 vlog_Y_a = -0.5*nvar*fss*log(2*pi) + fss*log(abs(det(A0xhat))) + Yexpt 0088 % a: given a0 and a+ 0089 logMarLHres = 0; % Initialize log of the marginal likelihood (restricted or constant parameters). 0090 for ki=1:fss %ndobs+1:fss % Forward recursion to get the marginal likelihood. See F on p.19 and pp. 48-49. 0091 %---- Restricted log marginal likelihood function (constant parameters). 0092 [A0l,A0u] = lu(A0xhat); 0093 ada = sum(log(abs(diag(A0u)))); % log|A0| 0094 termexp = y(ki,:)*A0xhat - phi(ki,:)*Apxhat; % 1-by-nvar 0095 logMarLHres = logMarLHres - (0.5*nvar)*log(2*pi) + ada - 0.5*termexp*termexp'; % log value 0096 end 0097 logMarLHres 0098 0099 0100 %--------- The log value of p(A+|Y,A0) at some point such as the peak ---------- 0101 totparsp = 0.0; 0102 tmpd = 0.0; 0103 for k=1:nvar 0104 totparsp = totparsp + size(Viconst{k},2); 0105 tmpd = tmpd + 0.5*log(abs(det(Hpinv{k}))); 0106 end 0107 vlog_ap_Ya0 = -0.5*totparsp*log(2*pi) + tmpd + Apexpt; 0108 0109 0110 0111 0112 %=================================== 0113 % Compute p(a0,k|Y,ao) at some point such as the peak (in this situation, we simply 0114 % generate results from the original Gibbs sampler). See FORECAST (2) pp.70-71 0115 %=================================== 0116 %--- Global set up for Gibbs. 0117 [Tinv,UT] = fn_gibbsrvar_setup(H0inv, Uiconst, Hpinv, Pmat, Viconst, nvar, fss); 0118 % 0119 vlog_a0_Yao = zeros(nvar,1); 0120 % 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 0121 vlog=zeros(ndraws2,1); 0122 for k=1:nvar 0123 bk = Uiconst{k}'*A0xhat(:,k); 0124 indx_ks=[k:nvar]; % the columns that exclude 1-(k-1)th columns 0125 A0gbs0 = A0hat; % starting at some point such as the peak 0126 nk = n0(k); 0127 0128 if k<nvar 0129 %--------- The 1st set of draws to be tossed away. ------------------ 0130 for draws = 1:ndraws1 0131 if ~mod(draws,nbuffer) 0132 disp(' ') 0133 disp(sprintf('The %dth column or equation in A0 with %d 1st tossed-away draws in Gibbs',k,draws)) 0134 end 0135 A0gbs1 = fn_gibbsrvar(A0gbs0,UT,nvar,fss,n0,indx_ks); 0136 A0gbs0=A0gbs1; % repeat the Gibbs sampling 0137 end 0138 0139 0140 %--------- The 2nd set of draws to be used. ------------------ 0141 for draws = 1:ndraws2 0142 if ~mod(draws,nbuffer) 0143 disp(' ') 0144 disp(sprintf('The %dth column or equation in A0 with %d usable draws in Gibbs',k,draws)) 0145 end 0146 [A0gbs1, Wcell] = fn_gibbsrvar(A0gbs0,UT,nvar,fss,n0,indx_ks); 0147 %------ See p.71, Forecast (II). 0148 %------ Computing p(a0_k|Y,a_others) at some point such as the peak along the dimensions of indx_ks. 0149 Vk = Tinv{k}\Wcell{k}; %V_k on p.71 of Forecast (II). 0150 gbeta = Vk\bk; % inv(V_k)*b_k on p.71 of Forecast (II) where alpha_k = b_k in our notation. 0151 [Vtq,Vtr]=qr(Vk',0); %To get inv(V_k)'*inv(V_k) in (*) on p.71 of Forecast (II). 0152 % 0153 vlog(draws) = 0.5*(fss+nk)*log(fss)-log(abs(det(Vk)))-0.5*(nk-1)*log(2*pi)-... 0154 0.5*(fss+1)*log(2)-gammaln(0.5*(fss+1))+fss*log(abs(gbeta(1)))-... 0155 0.5*fss*bk'*(Vtr\(Vtr'\bk)); 0156 0157 A0gbs0=A0gbs1; % repeat the Gibbs sampling 0158 end 0159 vlogm=max(vlog); 0160 qlog=vlog-vlogm; 0161 vlogxhat=vlogm-log(ndraws2)+log(sum(exp(qlog))); 0162 vlog_a0_Yao(k) = vlogxhat; 0163 % 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 0164 else 0165 disp(' ') 0166 disp(sprintf('The last(6th) column or equation in A0 with no Gibbs draws')) 0167 [A0gbs1, Wcell] = fn_gibbsrvar(A0gbs0,UT,nvar,fss,n0,indx_ks) 0168 %------ See p.71, Forecast (II). 0169 %------ Computing p(a0_k|Y,a_others) at some point such as the peak along the dimensions of indx_ks. 0170 Vk = Tinv{k}\Wcell{k}; %V_k on p.71 of Forecast (II). 0171 gbeta = Vk\bk; % inv(V_k)*b_k on p.71 of Forecast (II) where alpha_k = b_k in our notation. 0172 [Vtq,Vtr]=qr(Vk',0); %To get inv(V_k)'*inv(V_k) in (*) on p.71 of Forecast (II). 0173 % 0174 vloglast = 0.5*(fss+nk)*log(fss)-log(abs(det(Vk)))-0.5*(nk-1)*log(2*pi)-... 0175 0.5*(fss+1)*log(2)-gammaln(0.5*(fss+1))+fss*log(abs(gbeta(1)))-... 0176 0.5*fss*bk'*(Vtr\(Vtr'\bk)); 0177 vlog_a0_Yao(k) = vloglast; 0178 end 0179 end 0180 ndraws2 0181 0182 disp('Prior pdf -- log(p(a0hat, a+hat)):'); 0183 vlog_a0p 0184 disp('LH pdf -- log(p(Y|a0hat, a+hat)):'); 0185 vlog_Y_a 0186 disp('Posterior Kernal -- logp(ahat) + logp(Y|ahat):'); 0187 vlog_Y_a + vlog_a0p 0188 disp('Posterior pdf -- log(p(a0_i_hat|a0_other_hat, Y)):'); 0189 vlog_a0_Yao 0190 disp('Posterior pdf -- log(p(aphat|a0hat, Y)):'); 0191 vlog_ap_Ya0 0192 0193 %--------- The value of marginal density p(Y) ---------- 0194 disp(' '); 0195 disp(' '); 0196 disp('************ Marginal Likelihood of Y or Marginal Data Density: ************'); 0197 vlogY = vlog_a0p+vlog_Y_a-sum(vlog_a0_Yao)-vlog_ap_Ya0