diff --git a/matlab/DsgeVarLikelihood.m b/matlab/DsgeVarLikelihood.m
index 6fa15817d41d0e1028ddf93d8c6f0c84687ac40d..8c34d5741e6ccb900bbb8ec185b3395add8b8a91 100644
--- a/matlab/DsgeVarLikelihood.m
+++ b/matlab/DsgeVarLikelihood.m
@@ -1,8 +1,8 @@
-function [fval,cost_flag,ys,trend_coeff,info] = DsgeVarLikelihood(xparam1,gend)
+function [fval,cost_flag,ys,trend_coeff,info,PHI,SIGMAu,tmp2] = DsgeVarLikelihood(xparam1,gend)
% stephane.adjemian@ens.fr [06-17-2005]
global bayestopt_ exo_nbr dr_ estim_params_ Sigma_e_ options_ xparam1_test trend_coeff_
-global dsge_prior_weight targ_ xparam_
+global dsge_prior_weight targ_ xparam_ dr_
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
@@ -39,7 +39,7 @@ if strcmpi(targ_,'deep')
end
elseif strcmpi(targ_,'lambda')
indx = strmatch('dsge_prior_weight',estim_params_.param_names,'exact');
- if ~isempty(indx)
+ if ~isempty(indx)
xparam1 = xparam_;
xparam1(indx) = xparam_tmp;
end
@@ -160,6 +160,10 @@ end
% 3. theorretical moments (second order)
%------------------------------------------------------------------------------
tmp = lyapunov_symm(T,R*Q*R');
+
+tmpbis = R*Q*R';
+tmpbis = tmpbis(bayestopt_.mf,bayestopt_.mf);
+
NumberOfObservedVariables = size(options_.varobs,1);
NumberOfLags = options_.varlag;
k = NumberOfObservedVariables*NumberOfLags ;
@@ -197,17 +201,17 @@ assignin('base','GYX',GYX);
if ~isinf(dsge_prior_weight)
SIGMAu = dsge_prior_weight*gend*TheoreticalAutoCovarianceOfTheObservedVariables(:,:,1) + mYY ;
- tmp = dsge_prior_weight*gend*GYX + mYX;
- SIGMAu = SIGMAu - tmp*inv(dsge_prior_weight*gend*GXX+mXX)*tmp';
- SIGMAu = SIGMAu/(gend*(dsge_prior_weight+1));
-
+ tmp1 = dsge_prior_weight*gend*GYX + mYX;
+ tmp2 = inv(dsge_prior_weight*gend*GXX+mXX);
+ SIGMAu = SIGMAu - tmp1*tmp2*tmp1';
+ SIGMAu = SIGMAu / (gend*(dsge_prior_weight+1));
+ PHI = tmp2*tmp1';
prodlng1 = sum(gammaln(.5*((1+dsge_prior_weight)*gend- ...
NumberOfObservedVariables*NumberOfLags ...
+1-(1:NumberOfObservedVariables)')));
prodlng2 = sum(gammaln(.5*(dsge_prior_weight*gend- ...
NumberOfObservedVariables*NumberOfLags ...
- +1-(1:NumberOfObservedVariables)')));
-
+ +1-(1:NumberOfObservedVariables)')));
lik = .5*NumberOfObservedVariables*log(det(dsge_prior_weight*gend*GXX+mXX)) ...
+ .5*((dsge_prior_weight+1)*gend-k)*log(det((dsge_prior_weight+1)*gend*SIGMAu)) ...
- .5*NumberOfObservedVariables*log(det(dsge_prior_weight*gend*GXX)) ...
@@ -216,16 +220,14 @@ if ~isinf(dsge_prior_weight)
- .5*log(2)*NumberOfObservedVariables*((dsge_prior_weight+1)*gend-k) ...
+ .5*log(2)*NumberOfObservedVariables*(dsge_prior_weight*gend-k) ...
- prodlng1 + prodlng2;
-
-else
-
- tmp1 = inv(GXX)*GYX';
- tmp2 = GYY -GYX*tmp1;
-
- lik = -.5*sum(diag(inv(tmp2)*(mYY-2*tmp1'*mYX'+tmp1'*mXX*tmp1))) ...
+else % cod� par SM (s�rement pas exact... Que font ici les moments empiriques ?).
+ tmp1 = GYX;
+ tmp2 = inv(GXX);
+ PHI = tmp2*tmp1';
+ SIGMAu = GYY - tmp1*tmp2*tmp1;
+ % � finir de corriger...
+ lik = -.5*sum(diag(inv(tmp2)*(mYY-2*tmp1'*mYX'+tmp1'*mXX*tmp1))) ...
-(gend/2)*log(det(tmp2));
-
-
end
lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p1,bayestopt_.p2,bayestopt_.p3,bayestopt_.p4);
diff --git a/matlab/VarSampleMoments.m b/matlab/VarSampleMoments.m
index df8eb2a01d86ed799f58d00031ff7ad3b0266064..a76761e1024dd24dd4a0eeb39b2c639b686cef6a 100644
--- a/matlab/VarSampleMoments.m
+++ b/matlab/VarSampleMoments.m
@@ -49,7 +49,7 @@ for t=1:NumberOfObservations
elseif trend == 1
X(t,indx) = [ 1 , t ];
end
-end
+end
YtY = Y'*Y;
YtX = Y'*X;
diff --git a/matlab/dr1.m b/matlab/dr1.m
index 81dae7de1f8d984f4a76b26f40d51f09cf9ab229..520963d558db6f4411cf4b4f05da31c15b7936a5 100644
--- a/matlab/dr1.m
+++ b/matlab/dr1.m
@@ -620,11 +620,4 @@ end
% 06/01/2003 MJ added a test for xkmax_ > 1 and order > 1
% 08/28/2003 MJ corrected bug in computation of 2nd order (ordering of
% forward variable in LHS)
-% 08/29/2003 MJ use Sims routine if mjdgges isn't available
-
-
-
-
-
-
-
+% 08/29/2003 MJ use Sims routine if mjdgges isn't available
\ No newline at end of file
diff --git a/matlab/metropolis.m b/matlab/metropolis.m
index 11eb229f2b87d0f3e211667ee6c6d0598127cdb7..7e28c970f0a49a96e09d2bcf1fffdd3d8e3effa9 100644
--- a/matlab/metropolis.m
+++ b/matlab/metropolis.m
@@ -2370,7 +2370,6 @@ function metropolis(xparam1,vv,gend,data,rawdata,mh_bounds)
plot(1:10+options_.forecast,...
[data(idx,size(data,2)-10+1:end)';...
MeanForecast(:,k)],'-b','linewidth',2)
-
end
offsetx = 10;
else
@@ -2566,26 +2565,19 @@ function metropolis(xparam1,vv,gend,data,rawdata,mh_bounds)
end
% bvar-dsge
if ~isempty(dsge_prior_weight)
- DsgeVarLikelihood(deep',gend);
- GYY = evalin('base','GYY');
- GXX = evalin('base','GXX');
- GYX = evalin('base','GYX');
+ [fval,cost_flag,ys,trend_coeff,info,PHI,SIGMAu,iXX] = DsgeVarLikelihood(deep',gend);
DSGE_PRIOR_WEIGHT = floor(gend*(1+dsge_prior_weight));
- tmp1 = mYY + dsge_prior_weight*gend*GYY;
- tmp2 = mXY + dsge_prior_weight*gend*GYX';
- tmp3 = mXX + dsge_prior_weight*gend*GXX;
- tmp4 = (tmp1-tmp2'*inv(tmp3)*tmp2);
- PHI = inv(tmp3)*tmp2;
- val = 1;
+ tmp1 = SIGMAu*gend*(dsge_prior_weight+1);
+ val = 1;
+ tmp1 = chol(inv(tmp1))';
while val;
% draw from the marginal posterior of sig
- tmp5 = chol(inv(tmp4))'*randn(nvobs,DSGE_PRIOR_WEIGHT-NumberOfLagsTimesNvobs);
- SIGMAu_draw = inv(tmp5*tmp5');
- % draw from the conditional posterior of PHI
- VARvecPHI = kron(SIGMAu_draw,inv(tmp3));
+ tmp2 = tmp1*randn(nvobs,DSGE_PRIOR_WEIGHT-NumberOfLagsTimesNvobs);
+ SIGMAu_draw = inv(tmp2*tmp2');
+ % draw from the conditional posterior of PHI
+ VARvecPHI = kron(SIGMAu_draw,iXX);
PHI_draw = PHI(:) + chol(VARvecPHI)'*randn(nvobs*NumberOfLagsTimesNvobs,1);
- COMP_draw(1:nvobs,:) = reshape(PHI_draw, ...
- NumberOfLagsTimesNvobs,nvobs)';
+ COMP_draw(1:nvobs,:) = reshape(PHI_draw,NumberOfLagsTimesNvobs,nvobs)';
% Check for stationarity
tests = find(abs(eig(COMP_draw))>0.9999999999);
if isempty(tests)
@@ -2595,25 +2587,26 @@ function metropolis(xparam1,vv,gend,data,rawdata,mh_bounds)
% Get rotation
if dsge_prior_weight > 0
Atheta(dr_.order_var,:) = dr_.ghu*sqrt(Sigma_e_);
- A0 = Atheta(bayestopt_.mf,:);
+ A0 = Atheta(bayestopt_.mfys,:);
[OMEGAstar,SIGMAtr] = qr2(A0');
end
+ % SIGMAu_draw = A0*A0';
SIGMAu_chol = chol(SIGMAu_draw)';
SIGMAtrOMEGA = SIGMAu_chol*OMEGAstar';
PHIpower = eye(NumberOfLagsTimesNvobs);
irfs = zeros (nirfs,nvobs*exo_nbr);
- tmp6 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
- irfs(1,:) = tmp6(:)';
+ tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
+ irfs(1,:) = tmp3(:)';
for t = 2:nirfs
PHIpower = COMP_draw*PHIpower;
- tmp7 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
- irfs(t,:) = tmp7(:)';
+ tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
+ irfs(t,:) = tmp3(:)';
end
for j = 1:(nvobs*exo_nbr)
if max(irfs(:,j)) - min(irfs(:,j)) > 1e-10
tmp_dsgevar(:,j) = (irfs(:,j));
- end
- end
+ end
+ end
if irun_irf_dsgevar < MAX_nirfs_dsgevar
stock_irf_dsgevar(:,:,:,irun_irf_dsgevar) = reshape(tmp_dsgevar,nirfs,nvobs,exo_nbr);
else
@@ -2649,7 +2642,7 @@ function metropolis(xparam1,vv,gend,data,rawdata,mh_bounds)
irun_irf_dsge = 0;
if ~isempty(dsge_prior_weight)
sfil_irf_dsgevar = 0;
- irun_irf_dsgevar = 0;
+ irun_irf_dsgevar = 0;
end
eval(['load ' instr1 int2str(ffil) instr2]);
NumberOfSimulations = length(logpo2);
@@ -2695,26 +2688,19 @@ function metropolis(xparam1,vv,gend,data,rawdata,mh_bounds)
end
% bvar-dsge
if ~isempty(dsge_prior_weight)
- DsgeVarLikelihood(deep',gend);
- GYY = evalin('base','GYY');
- GXX = evalin('base','GXX');
- GYX = evalin('base','GYX');
+ [fval,cost_flag,ys,trend_coeff,info,PHI,SIGMAu,iXX] = DsgeVarLikelihood(deep',gend);
DSGE_PRIOR_WEIGHT = floor(gend*(1+dsge_prior_weight));
- tmp1 = mYY + dsge_prior_weight*gend*GYY;
- tmp2 = mXY + dsge_prior_weight*gend*GYX';
- tmp3 = mXX + dsge_prior_weight*gend*GXX;
- tmp4 = (tmp1-tmp2'*inv(tmp3)*tmp2);
- PHI = inv(tmp3)*tmp2;
- val = 1;
+ tmp1 = SIGMAu*gend*(1+dsge_prior_weight);
+ tmp1 = chol(inv(tmp1))';
+ val = 1;
while val;
% draw from the marginal posterior of sig
- tmp5 = chol(inv(tmp4))'*randn(nvobs,DSGE_PRIOR_WEIGHT-NumberOfLagsTimesNvobs);
- SIGMAu_draw = inv(tmp5*tmp5');
+ tmp2 = tmp1*randn(nvobs,DSGE_PRIOR_WEIGHT-NumberOfLagsTimesNvobs);
+ SIGMAu_draw = inv(tmp2*tmp2');
% draw from the conditional posterior of PHI
- VARvecPHI = kron(SIGMAu_draw,inv(tmp3));
+ VARvecPHI = kron(SIGMAu_draw,iXX);
PHI_draw = PHI(:) + chol(VARvecPHI)'*randn(nvobs*NumberOfLagsTimesNvobs,1);
- COMP_draw(1:nvobs,:) = reshape(PHI_draw, ...
- NumberOfLagsTimesNvobs,nvobs)';
+ COMP_draw(1:nvobs,:) = reshape(PHI_draw,NumberOfLagsTimesNvobs,nvobs)';
% Check for stationarity
tests = find(abs(eig(COMP_draw))>0.9999999999);
if isempty(tests)
@@ -2724,20 +2710,20 @@ function metropolis(xparam1,vv,gend,data,rawdata,mh_bounds)
% Get rotation
if dsge_prior_weight > 0
Atheta(dr_.order_var,:) = dr_.ghu*sqrt(Sigma_e_);
- A0 = Atheta(bayestopt_.mf,:);
+ A0 = Atheta(bayestopt_.mfys,:);
[OMEGAstar,SIGMAtr] = qr2(A0');
end
SIGMAu_chol = chol(SIGMAu_draw)';
SIGMAtrOMEGA = SIGMAu_chol*OMEGAstar';
PHIpower = eye(NumberOfLagsTimesNvobs);
irfs = zeros (nirfs,nvobs*exo_nbr);
- tmp6 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
- irfs(1,:) = tmp6(:)';
+ tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
+ irfs(1,:) = tmp3(:)';
for t = 2:nirfs
PHIpower = COMP_draw*PHIpower;
- tmp7 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
- irfs(t,:) = tmp7(:)';
- end
+ tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
+ irfs(t,:) = tmp3(:)';
+ end
for j = 1:(nvobs*exo_nbr)
if max(irfs(:,j)) - min(irfs(:,j)) > 1e-10
tmp_dsgevar(:,j) = (irfs(:,j));
diff --git a/matlab/metropolis99.m b/matlab/metropolis99.m
index 53670b0c4dbdf70fa267b25295b831f575304022..5d8edbbc0a696520599d6d2140a3ce32b89e43f5 100644
--- a/matlab/metropolis99.m
+++ b/matlab/metropolis99.m
@@ -1,17 +1,15 @@
function [PostMod,PostVar,Scale,PostMean] = ...
samaxwell(xparam1,gend,data,mh_bounds,NumberOfIterationsForInitialization,init_scale,info,MeanPar,VarCov)
% stephane.adjemian@cepremap.cnrs.fr [11-22-2005]
-% Adapted from an older version of metropolis.m
-
global bayestopt_ exo_nbr dr_ estim_params_ Sigma_e_ options_ xparam_test
global lgy_ lgx_ fname_ ys_ xkmin_ xkmax_ ykmin_ ykmax_ endo_nbr mean_varobs
global oo_ lgx_orig_ord_ lgy_TeX_ lgx_TeX_ dsge_prior_weight
-
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
ncn = estim_params_.ncn;
np = estim_params_.np ;
+ns = nvx+nvn+ncx+ncn;
nx = nvx+nvn+ncx+ncn+np;
npar = length(xparam1);
nvobs = size(options_.varobs,1);
@@ -96,7 +94,7 @@ set(hh,'Name','Looking for the posterior covariance...')
j = 1;
isux = 0;
% ix2 = xparam1;% initial condition!
-if isempty(strmatch('dsge_prior_weight',estim_params_.param_names)) & isempty(dsge_prior_weight)
+if isempty(strmatch('dsge_prior_weight',estim_params_.param_names)) & Isempty(dsge_prior_weight)
ilogpo2 = -DsgeLikelihood(ix2,gend,data);
else
ilogpo2 = -DsgeVarLikelihood(ix2,gend);
@@ -178,7 +176,7 @@ if strcmpi(info,'LastCall')
jsux = jsux + 1;
else
% ... otherwise I don't move.
- end
+ end
prtfrc = j/MaxNumberOfTuningSimulations;
waitbar(prtfrc,hh,sprintf('%f done, acceptation rate %f',prtfrc,isux/j));
if j/500 == round(j/500)
@@ -204,18 +202,14 @@ if strcmpi(info,'LastCall')
%%
%% Now I climb the hill
%%
- hh = waitbar(0,'');
+ hh = waitbar(0,' ');
set(hh,'Name','Now I am climbing the hill...')
j = 1; jj = 1;
jsux = 0;
test = 0;
- % ix2 = ModePar;%% I start from the point with the highest posterior density...
- % if isempty(strmatch('dsge_prior_weight',estim_params_.param_names)) & isempty(dsge_prior_weight)
- % ilogpo2 = -DsgeLikelihood(ix2,gend,data);% initial posterior density
- % else
- % ilogpo2 = -DsgeVarLikelihood(ix2,gend);% initial posterior density
- % end
dd = transpose(chol(CovJump));
+ init_scale = eye(length(dd))*init_scale;
+ indx = strmatch('dsge_prior_weight',estim_params_.param_names,'exact');
while j<=MaxNumberOfClimbingSimulations
proposal = init_scale*dd*randn(npar,1) + ModePar;
if all(proposal > mh_bounds(:,1)) & all(proposal < mh_bounds(:,2))
@@ -242,13 +236,18 @@ if strcmpi(info,'LastCall')
%init_scale = init_scale*cfactor;
jsux = 0;
jj = 0;
- if test1<0.005
+ if test1<0.001
test = test+1;
end
if test>4% If I do not progress enough I reduce the scale parameter
% of the jumping distribution (cooling down the system).
-
- init_scale = init_scale/1.20;
+ if isempty(indx)
+ init_scale = init_scale/1.10;
+ else
+ init_scale = init_scale/1.10;
+ % USEFUL IF THE MODEL IS TRUE:
+ % init_scale(indx+ns,indx+ns) = init_scale(indx+ns,indx+ns)*1.10/1.001;
+ end
end
end
j = j+1;
diff --git a/matlab/subset.m b/matlab/subset.m
index 41eae520a1c8e951b7b8d37359476dd87f9b1bb3..a82ab1c9774a574ca8f2e0c088b26c24bfb06cfc 100644
--- a/matlab/subset.m
+++ b/matlab/subset.m
@@ -7,7 +7,6 @@ VarObs = options_.varobs;
VarExo = lgx_;
info = options_.ParamSubSet;
-
nvx = estim_params_.nvx;
nvn = estim_params_.nvn;
ncx = estim_params_.ncx;
@@ -27,8 +26,10 @@ elseif strcmpi(info,'MeasurementErrors')
indx = [(nvx+1:nvx+nvn),(nvx+nvn+ncx+1:nvx+nvn+ncx+ncn)]';
elseif strcmpi(info,'MeasurementErrorsWithoutCorrelations')
indx = (nvx+1:nvx+nvn)';
-elseif strcmpi(info,'AllWithoutMeasurementErros')
+elseif strcmpi(info,'AllWithoutMeasurementErrors')
indx = [(1:nvx),nvx+nvn+1:nvx+nvn+ncx,nvx+nvn+ncx+ncn+1:nx]';
+elseif strcmpi(info,'None')
+ indx = [];
end
if isempty(ExcludedParamNames)