stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Factorization of rand_inverse_wishart(...) and rand_matrix_normal(...) + Some small changes to (marginally) speed up the thing (loop --> matrix calculus in rand_inverse_wishart(...), ... ). 

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1398 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
adjemian 2007-09-25 16:25:25 +00:00
parent 0c3d36eb11
commit 4a0aeb12de
4 changed files with 76 additions and 81 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
matlab/PosteriorIRF.m
View File

@ -123,7 +123,7 @@ if MAX_nirfs_dsgevar
var_sample_moments(options_.first_obs,options_.first_obs+options_.nobs-1,options_.varlag,-1);
NumberOfLags = options_.varlag;
NumberOfLagsTimesNvobs = NumberOfLags*nvobs;
COMP_draw = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs);
Companion_matrix = diag(ones(nvobs*(NumberOfLags-1),1),-nvobs);
end
b = 0;
nosaddle = 0;
@ -179,22 +179,18 @@ while b<=B
[fval,cost_flag,ys,trend_coeff,info,PHI,SIGMAu,iXX] = DsgeVarLikelihood(deep',gend);
dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names));
DSGE_PRIOR_WEIGHT = floor(gend*(1+dsge_prior_weight));
tmp1 = SIGMAu*gend*(dsge_prior_weight+1);
val = 1;
tmp1 = chol(inv(tmp1))';
while val;
% draw from the marginal posterior of sig
tmp2 = tmp1*randn(nvobs,DSGE_PRIOR_WEIGHT-NumberOfLagsTimesNvobs);
SIGMAu_draw = inv(tmp2*tmp2');
SIGMA_inv_upper_chol = chol(inv(SIGMAu*gend*(dsge_prior_weight+1)));
explosive_var = 1;
while explosive_var
% draw from the marginal posterior of SIGMA
SIGMAu_draw = rand_inverse_wishart(nvobs, DSGE_PRIOR_WEIGHT-NumberOfLagsTimesNvobs, ...
SIGMA_inv_upper_chol);
% 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)';
PHI_draw = rand_matrix_normal(NumberOfLagsTimesNvobs,nvobs, PHI, ...
chol(iXX)', chol(SIGMAu_draw)');
Companion_matrix(1:nvobs,:) = transpose(PHI_draw);
% Check for stationarity
tests = find(abs(eig(COMP_draw))>0.9999999999);
if isempty(tests)
val=0;
end
explosive_var = any(abs(eig(Companion_matrix))>1.000000001);
end
% Get rotation
if dsge_prior_weight > 0
@ -209,7 +205,7 @@ while b<=B
tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
irfs(1,:) = tmp3(:)';
for t = 2:options_.irf
PHIpower = COMP_draw*PHIpower;
PHIpower = Companion_matrix*PHIpower;
tmp3 = PHIpower(1:nvobs,1:nvobs)*SIGMAtrOMEGA;
irfs(t,:) = tmp3(:)';
end
@ -222,7 +218,8 @@ while b<=B
stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,nvobs,M_.exo_nbr);
else
stock_irf_bvardsge(:,:,:,IRUN) = reshape(tmp_dsgevar,options_.irf,nvobs,M_.exo_nbr);
instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) ' stock_irf_bvardsge;'];,
instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' ...
int2str(NumberOfIRFfiles_dsgevar) ' stock_irf_bvardsge;'];,
eval(['save ' instr]);
NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
IRUN =0;
@ -234,7 +231,8 @@ while b<=B
stock_irf_dsge = stock_irf_dsge(:,:,:,1:irun);
if MAX_nirfs_dsgevar & (b == B | IRUN == B)
stock_irf_bvardsge = stock_irf_bvardsge(:,:,:,1:IRUN);
instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' int2str(NumberOfIRFfiles_dsgevar) ' stock_irf_bvardsge;'];,
instr = [MhDirectoryName '/' M_.fname '_irf_bvardsge' ...
int2str(NumberOfIRFfiles_dsgevar) ' stock_irf_bvardsge;'];,
eval(['save ' instr]);
NumberOfIRFfiles_dsgevar = NumberOfIRFfiles_dsgevar+1;
irun = 0;

71
matlab/bvar_forecast.m
View File

@ -46,10 +46,12 @@ function bvar_forecast(nlags)
% Without shocks
lags_data = forecast_data.initval;
for t = 1:options_.forecast
X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.xdata(t, :) ];
y = X * Phi;
lags_data = [ lags_data(2:end, :); y ];
lags_data(1:end-1,:) = lags_data(2:end, :);
lags_data(end,:) = y;
sims_no_shock(t, :, d) = y;
end
@ -59,7 +61,8 @@ function bvar_forecast(nlags)
X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.xdata(t, :) ];
shock = (Sigma_lower_chol * randn(ny, 1))';
y = X * Phi + shock;
lags_data = [ lags_data(2:end, :); y ];
lags_data(1:end-1,:) = lags_data(2:end, :);
lags_data(end,:) = y;
sims_with_shocks(t, :, d) = y;
end
end
@ -104,7 +107,8 @@ function bvar_forecast(nlags)
for t = 1:size(forecast_data.realized_val, 1)
X = [ reshape(flipdim(lags_data, 1)', 1, ny*nlags) forecast_data.realized_xdata(t, :) ];
y = X * posterior.PhiHat;
lags_data = [ lags_data(2:end, :); y ];
lags_data(1:end-1,:) = lags_data(2:end, :);
lags_data(end,:) = y;
sq_err_cumul = sq_err_cumul + (y - forecast_data.realized_val(t, :)) .^ 2;
end
@ -115,63 +119,4 @@ function bvar_forecast(nlags)
for i = 1:size(options_.varobs, 1)
fprintf('%s: %10.4f\n', options_.varobs(i, :), rmse(i));
end
end
function G = rand_inverse_wishart(m, v, H_inv_upper_chol)
% rand_inverse_wishart Pseudo random matrices drawn from an
% inverse Wishart distribution
%
% G = rand_inverse_wishart(m, v, H_inv_upper_chol)
%
% Returns an m-by-m matrix drawn from an inverse-Wishart distribution.
%
% m: dimension of G and H_inv_upper_chol.
% v: degrees of freedom, greater or equal than m.
% H_inv_chol: upper cholesky decomposition of the inverse of the
% matrix parameter.
% The upper cholesky of the inverse is requested here
% in order to avoid to recompute it at every random draw.
% H_inv_upper_chol = chol(inv(H))
%
% In other words:
% G ~ IW(m, v, H) where H = inv(H_inv_upper_chol'*H_inv_upper_chol)
% or, equivalently, using the correspondence between Wishart and
% inverse-Wishart:
% inv(G) ~ W(m, v, S) where S = H_inv_upper_chol'*H_inv_upper_chol = inv(H)
X = NaN(v, m);
for i = 1:v
X(i, :) = randn(1, m) * H_inv_upper_chol;
end
% At this point, X'*X is Wishart distributed
% G = inv(X'*X);
% Rather compute inv(X'*X) using the SVD
[U,S,V] = svd(X, 0);
SSi = 1 ./ (diag(S) .^ 2);
G = (V .* repmat(SSi', m, 1)) * V';
function B = rand_matrix_normal(n, p, M, Omega_lower_chol, Sigma_lower_chol)
% rand_matrix_normal Pseudo random matrices drawn from a
% matrix-normal distribution
%
% B = rand_matrix_normal(n, p, M, Omega_lower_chol, Sigma_lower_chol)
%
% Returns an n-by-p matrix drawn from a Matrix-normal distribution
% Same notations than: http://en.wikipedia.org/wiki/Matrix_normal_distribution
% M is the mean, n-by-p matrix
% Omega_lower_chol is n-by-n, lower Cholesky decomposition of Omega
% (Omega_lower_chol = chol(Omega, 'lower'))
% Sigma_lower_chol is p-by-p, lower Cholesky decomposition of Sigma
% (Sigma_lower_chol = chol(Sigma, 'lower'))
%
% Equivalent to vec(B) ~ N(vec(Mu), kron(Sigma, Omega))
%
B1 = randn(n * p, 1);
B2 = kron(Sigma_lower_chol, Omega_lower_chol) * B1;
B3 = reshape(B2, n, p);
B = B3 + M;
end

32
matlab/rand_inverse_wishart.m Normal file
View File

@ -0,0 +1,32 @@
function G = rand_inverse_wishart(m, v, H_inv_upper_chol)
% rand_inverse_wishart Pseudo random matrices drawn from an
% inverse Wishart distribution
%
% G = rand_inverse_wishart(m, v, H_inv_upper_chol)
%
% Returns an m-by-m matrix drawn from an inverse-Wishart distribution.
%
% m: dimension of G and H_inv_upper_chol.
% v: degrees of freedom, greater or equal than m.
% H_inv_chol: upper cholesky decomposition of the inverse of the
% matrix parameter.
% The upper cholesky of the inverse is requested here
% in order to avoid to recompute it at every random draw.
% H_inv_upper_chol = chol(inv(H))
%
% In other words:
% G ~ IW(m, v, H) where H = inv(H_inv_upper_chol'*H_inv_upper_chol)
% or, equivalently, using the correspondence between Wishart and
% inverse-Wishart:
% inv(G) ~ W(m, v, S) where S = H_inv_upper_chol'*H_inv_upper_chol = inv(H)
X = randn(v, m) * H_inv_upper_chol;
% At this point, X'*X is Wishart distributed
% G = inv(X'*X);
% Rather compute inv(X'*X) using the SVD
[U,S,V] = svd(X, 0);
SSi = 1 ./ (diag(S) .^ 2);
G = (V .* repmat(SSi', m, 1)) * V';

20
matlab/rand_matrix_normal.m Normal file
View File

@ -0,0 +1,20 @@
function B = rand_matrix_normal(n, p, M, Omega_lower_chol, Sigma_lower_chol)
% rand_matrix_normal Pseudo random matrices drawn from a
% matrix-normal distribution
%
% B = rand_matrix_normal(n, p, M, Omega_lower_chol, Sigma_lower_chol)
%
% Returns an n-by-p matrix drawn from a Matrix-normal distribution
% Same notations than: http://en.wikipedia.org/wiki/Matrix_normal_distribution
% M is the mean, n-by-p matrix
% Omega_lower_chol is n-by-n, lower Cholesky decomposition of Omega
% (Omega_lower_chol = chol(Omega, 'lower'))
% Sigma_lower_chol is p-by-p, lower Cholesky decomposition of Sigma
% (Sigma_lower_chol = chol(Sigma, 'lower'))
%
% Equivalent to vec(B) ~ N(vec(Mu), kron(Sigma, Omega))
B1 = randn(n * p, 1);
B2 = kron(Sigma_lower_chol, Omega_lower_chol) * B1;
B3 = reshape(B2, n, p);
B = B3 + M;