Estimation C++ DLL: 4 fixes and some additional files used in testing: sweuromodel_dll.mod, data and DsgeLikelihood.m modified for testing under octave

time-shift
George Perendia 2010-07-23 12:08:00 +01:00
parent a0dc2ce4ee
commit 1ac71f3d0c
5 changed files with 1507 additions and 6 deletions

View File

@ -65,7 +65,7 @@ KalmanFilter::compute_zeta_varobs_back_mixed(const std::vector<size_t> &zeta_bac
zeta_mixed_arg.begin(), zeta_mixed_arg.end(),
back_inserter(zeta_back_mixed));
set_union(zeta_back_mixed.begin(), zeta_back_mixed.end(),
varobs_arg.begin(), varobs_arg.end(),
varobs_sorted.begin(), varobs_sorted.end(),
back_inserter(zeta_varobs_back_mixed));
return zeta_varobs_back_mixed;
}
@ -95,7 +95,6 @@ KalmanFilter::filter(const MatrixView &detrendedDataView, const Matrix &H, Vect
{
double loglik=0.0, ll, logFdet, Fdet;
size_t p = Finv.getRows();
bool nonstationary = true;
for (size_t t = 0; t < detrendedDataView.getCols(); ++t)
{
@ -130,7 +129,12 @@ KalmanFilter::filter(const MatrixView &detrendedDataView, const Matrix &H, Vect
blas::gemm("N", "N", 1.0, T, Pstar, 0.0, Ptmp);
// 3) Pt+1= Ptmp*T' +RQR'
Pstar = RQRt;
blas::gemm("N", "T", 1.0, Ptmp, T, 1.0, Pstar);
//blas::gemm("N", "T", 1.0, Ptmp, T, 1.0, Pstar);
//enforce Pstar symmetry with P=(P+P')/2=0.5P+0.5P'
blas::gemm("N", "T", 0.5, Ptmp, T, 0.5, Pstar);
mat::transpose(Ptmp, Pstar);
mat::add(Pstar,Ptmp);
if (t>0)
nonstationary = mat::isDiff(KFinv, oldKFinv, riccati_tol);
oldKFinv=KFinv;
@ -139,6 +143,7 @@ KalmanFilter::filter(const MatrixView &detrendedDataView, const Matrix &H, Vect
// err= Yt - Za
MatrixConstView yt(detrendedDataView, 0, t, detrendedDataView.getRows(), 1); // current observation vector
vt = yt;
blas::gemm("N", "N", -1.0, Z, a_init, 1.0, vt);
// at+1= T(at+ KFinv *err)
blas::gemm("N", "N", 1.0, KFinv, vt, 1.0, a_init);

View File

@ -132,6 +132,7 @@ logposterior(const VectorConstView &estParams, const MatrixConstView &data,
double qz_criterium = *mxGetPr(mxGetField(options_, 0, "qz_criterium"));
double lyapunov_tol = *mxGetPr(mxGetField(options_, 0, "lyapunov_complex_threshold"));
double riccati_tol = *mxGetPr(mxGetField(options_, 0, "riccati_tol"));
size_t presample = (size_t) *mxGetPr(mxGetField(options_, 0, "presample"));
std::vector<size_t> varobs;
const mxArray *varobs_mx = mxGetField(options_, 0, "varobs_id");
@ -175,16 +176,17 @@ logposterior(const VectorConstView &estParams, const MatrixConstView &data,
Vector deepParams(n_param);
deepParams = VectorConstView(mxGetPr(mxGetField(M_, 0, "params")), n_param, 1);
Matrix Q(n_exo);
Q = MatrixConstView(mxGetPr(mxGetField(M_, 0, "Sigma_e")), n_exo, n_exo, 1);
Q = MatrixConstView(mxGetPr(mxGetField(M_, 0, "Sigma_e")), n_exo, n_exo, n_exo);
Matrix H(n_varobs);
const mxArray *H_mx = mxGetField(M_, 0, "H");
if (mxGetM(H_mx) == 1 && mxGetN(H_mx) == 1 && *mxGetPr(H_mx) == 0)
H.setAll(0.0);
else
H = MatrixConstView(mxGetPr(mxGetField(M_, 0, "H")), n_varobs, n_varobs, 1);
H = MatrixConstView(mxGetPr(mxGetField(M_, 0, "H")), n_varobs, n_varobs, n_varobs);
// Compute the posterior
double logPD = lpd.compute(steadyState, estParams2, deepParams, data, Q, H, 0, info);
double logPD =lpd.compute(steadyState, estParams2, deepParams, data, Q, H, presample, info);
// Cleanups
for (std::vector<EstimatedParameter>::iterator it = estParamsInfo.begin();

View File

@ -0,0 +1,345 @@
function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
% function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)
% Evaluates the posterior kernel of a dsge model.
%
% INPUTS
% xparam1 [double] vector of model parameters.
% gend [integer] scalar specifying the number of observations.
% data [double] matrix of data
% data_index [cell] cell of column vectors
% number_of_observations [integer]
% no_more_missing_observations [integer]
% OUTPUTS
% fval : value of the posterior kernel at xparam1.
% cost_flag : zero if the function returns a penalty, one otherwise.
% ys : steady state of original endogenous variables
% trend_coeff :
% info : vector of informations about the penalty:
% 41: one (many) parameter(s) do(es) not satisfied the lower bound
% 42: one (many) parameter(s) do(es) not satisfied the upper bound
%
% SPECIAL REQUIREMENTS
%
% Copyright (C) 2004-2009, 2010 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/>.
global bayestopt_ estim_params_ options_ trend_coeff_ M_ oo_
fval = [];
ys = [];
trend_coeff = [];
cost_flag = 1;
nobs = size(options_.varobs,1);
%------------------------------------------------------------------------------
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
if options_.mode_compute ~= 1 & any(xparam1 < bayestopt_.lb)
k = find(xparam1 < bayestopt_.lb);
fval = bayestopt_.penalty+sum((bayestopt_.lb(k)-xparam1(k)).^2);
cost_flag = 0;
info = 41;
return;
end
if options_.mode_compute ~= 1 & any(xparam1 > bayestopt_.ub)
k = find(xparam1 > bayestopt_.ub);
fval = bayestopt_.penalty+sum((xparam1(k)-bayestopt_.ub(k)).^2);
cost_flag = 0;
info = 42;
return;
end
Q = M_.Sigma_e;
H = M_.H;
for i=1:estim_params_.nvx
k =estim_params_.var_exo(i,1);
Q(k,k) = xparam1(i)*xparam1(i);
end
offset = estim_params_.nvx;
if estim_params_.nvn
for i=1:estim_params_.nvn
k = estim_params_.var_endo(i,1);
H(k,k) = xparam1(i+offset)*xparam1(i+offset);
end
offset = offset+estim_params_.nvn;
end
if estim_params_.ncx
for i=1:estim_params_.ncx
k1 =estim_params_.corrx(i,1);
k2 =estim_params_.corrx(i,2);
Q(k1,k2) = xparam1(i+offset)*sqrt(Q(k1,k1)*Q(k2,k2));
Q(k2,k1) = Q(k1,k2);
end
[CholQ,testQ] = chol(Q);
if testQ %% The variance-covariance matrix of the structural innovations is not definite positive.
%% We have to compute the eigenvalues of this matrix in order to build the penalty.
a = diag(eig(Q));
k = find(a < 0);
if k > 0
fval = bayestopt_.penalty+sum(-a(k));
cost_flag = 0;
info = 43;
return
end
end
offset = offset+estim_params_.ncx;
end
if estim_params_.ncn
for i=1:estim_params_.ncn
k1 = options_.lgyidx2varobs(estim_params_.corrn(i,1));
k2 = options_.lgyidx2varobs(estim_params_.corrn(i,2));
H(k1,k2) = xparam1(i+offset)*sqrt(H(k1,k1)*H(k2,k2));
H(k2,k1) = H(k1,k2);
end
[CholH,testH] = chol(H);
if testH
a = diag(eig(H));
k = find(a < 0);
if k > 0
fval = bayestopt_.penalty+sum(-a(k));
cost_flag = 0;
info = 44;
return
end
end
offset = offset+estim_params_.ncn;
end
if estim_params_.np > 0
M_.params(estim_params_.param_vals(:,1)) = xparam1(offset+1:end);
end
M_.Sigma_e = Q;
M_.H = H;
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
[T,R,SteadyState,info] = dynare_resolve(bayestopt_.restrict_var_list,...
bayestopt_.restrict_columns,...
bayestopt_.restrict_aux);
if info(1) == 1 || info(1) == 2 || info(1) == 5
fval = bayestopt_.penalty+1;
cost_flag = 0;
return
elseif info(1) == 3 || info(1) == 4 || info(1)==6 ||info(1) == 19 || info(1) == 20 || info(1) == 21
fval = bayestopt_.penalty+info(2);
cost_flag = 0;
return
end
bayestopt_.mf = bayestopt_.mf1;
if options_.noconstant
constant = zeros(nobs,1);
else
if options_.loglinear
constant = log(SteadyState(bayestopt_.mfys));
else
constant = SteadyState(bayestopt_.mfys);
end
end
if bayestopt_.with_trend
trend_coeff = zeros(nobs,1);
t = options_.trend_coeffs;
for i=1:length(t)
if ~isempty(t{i})
trend_coeff(i) = evalin('base',t{i});
end
end
trend = repmat(constant,1,gend)+trend_coeff*[1:gend];
else
trend = repmat(constant,1,gend);
end
start = options_.presample+1;
np = size(T,1);
mf = bayestopt_.mf;
no_missing_data_flag = (number_of_observations==gend*nobs);
%------------------------------------------------------------------------------
% 3. Initial condition of the Kalman filter
%------------------------------------------------------------------------------
T
R
Q
R*Q*R'
pause
options_.lik_init = 1;
kalman_algo = options_.kalman_algo;
if options_.lik_init == 1 % Kalman filter
if kalman_algo ~= 2
kalman_algo = 1;
end
Pstar = lyapunov_symm(T,R*Q*R',options_.qz_criterium,options_.lyapunov_complex_threshold);
Pinf = [];
elseif options_.lik_init == 2 % Old Diffuse Kalman filter
if kalman_algo ~= 2
kalman_algo = 1;
end
Pstar = options_.Harvey_scale_factor*eye(np);
Pinf = [];
elseif options_.lik_init == 3 % Diffuse Kalman filter
if kalman_algo ~= 4
kalman_algo = 3;
end
[QT,ST] = schur(T);
e1 = abs(ordeig(ST)) > 2-options_.qz_criterium;
[QT,ST] = ordschur(QT,ST,e1);
k = find(abs(ordeig(ST)) > 2-options_.qz_criterium);
nk = length(k);
nk1 = nk+1;
Pinf = zeros(np,np);
Pinf(1:nk,1:nk) = eye(nk);
Pstar = zeros(np,np);
B = QT'*R*Q*R'*QT;
for i=np:-1:nk+2
if ST(i,i-1) == 0
if i == np
c = zeros(np-nk,1);
else
c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
end
q = eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i);
Pstar(nk1:i,i) = q\(B(nk1:i,i)+c);
Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
else
if i == np
c = zeros(np-nk,1);
c1 = zeros(np-nk,1);
else
c = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i,i+1:end)')+...
ST(i,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i)+...
ST(i,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1);
c1 = ST(nk1:i,:)*(Pstar(:,i+1:end)*ST(i-1,i+1:end)')+...
ST(i-1,i-1)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i-1)+...
ST(i-1,i)*ST(nk1:i,i+1:end)*Pstar(i+1:end,i);
end
q = [eye(i-nk)-ST(nk1:i,nk1:i)*ST(i,i) -ST(nk1:i,nk1:i)*ST(i,i-1);...
-ST(nk1:i,nk1:i)*ST(i-1,i) eye(i-nk)-ST(nk1:i,nk1:i)*ST(i-1,i-1)];
z = q\[B(nk1:i,i)+c;B(nk1:i,i-1)+c1];
Pstar(nk1:i,i) = z(1:(i-nk));
Pstar(nk1:i,i-1) = z(i-nk+1:end);
Pstar(i,nk1:i-1) = Pstar(nk1:i-1,i)';
Pstar(i-1,nk1:i-2) = Pstar(nk1:i-2,i-1)';
i = i - 1;
end
end
if i == nk+2
c = ST(nk+1,:)*(Pstar(:,nk+2:end)*ST(nk1,nk+2:end)')+ST(nk1,nk1)*ST(nk1,nk+2:end)*Pstar(nk+2:end,nk1);
Pstar(nk1,nk1)=(B(nk1,nk1)+c)/(1-ST(nk1,nk1)*ST(nk1,nk1));
end
Z = QT(mf,:);
R1 = QT'*R;
[QQ,RR,EE] = qr(Z*ST(:,1:nk),0);
k = find(abs(diag([RR; zeros(nk-size(Z,1),size(RR,2))])) < 1e-8);
if length(k) > 0
k1 = EE(:,k);
dd =ones(nk,1);
dd(k1) = zeros(length(k1),1);
Pinf(1:nk,1:nk) = diag(dd);
end
end
if kalman_algo == 2
end
kalman_tol = options_.kalman_tol;
riccati_tol = options_.riccati_tol;
mf = bayestopt_.mf1;
Y = data-trend;
Pstar
pause
%------------------------------------------------------------------------------
% 4. Likelihood evaluation
%------------------------------------------------------------------------------
if (kalman_algo==1)% Multivariate Kalman Filter
if no_missing_data_flag
LIK = kalman_filter(T,R,Q,H,Pstar,Y,start,mf,kalman_tol,riccati_tol);
else
LIK = ...
missing_observations_kalman_filter(T,R,Q,H,Pstar,Y,start,mf,kalman_tol,riccati_tol, ...
data_index,number_of_observations,no_more_missing_observations);
end
if isinf(LIK)
kalman_algo = 2;
end
end
if (kalman_algo==2)% Univariate Kalman Filter
no_correlation_flag = 1;
if length(H)==1 & H == 0
H = zeros(nobs,1);
else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
H = diag(H);
else
no_correlation_flag = 0;
end
end
if no_correlation_flag
LIK = univariate_kalman_filter(T,R,Q,H,Pstar,Y,start,mf,kalman_tol,riccati_tol,data_index,number_of_observations,no_more_missing_observations);
else
LIK = univariate_kalman_filter_corr(T,R,Q,H,Pstar,Y,start,mf,kalman_tol,riccati_tol,data_index,number_of_observations,no_more_missing_observations);
end
end
if (kalman_algo==3)% Multivariate Diffuse Kalman Filter
if no_missing_data_flag
LIK = diffuse_kalman_filter(ST,R1,Q,H,Pinf,Pstar,Y,start,Z,kalman_tol, ...
riccati_tol);
else
LIK = missing_observations_diffuse_kalman_filter(ST,R1,Q,H,Pinf, ...
Pstar,Y,start,Z,kalman_tol,riccati_tol,...
data_index,number_of_observations,...
no_more_missing_observations);
end
if isinf(LIK)
kalman_algo = 4;
end
end
if (kalman_algo==4)% Univariate Diffuse Kalman Filter
no_correlation_flag = 1;
if length(H)==1 & H == 0
H = zeros(nobs,1);
else
if all(all(abs(H-diag(diag(H)))<1e-14))% ie, the covariance matrix is diagonal...
H = diag(H);
else
no_correlation_flag = 0;
end
end
if no_correlation_flag
LIK = univariate_diffuse_kalman_filter(ST,R1,Q,H,Pinf,Pstar,Y, ...
start,Z,kalman_tol,riccati_tol,data_index,...
number_of_observations,no_more_missing_observations);
else
LIK = univariate_diffuse_kalman_filter_corr(ST,R1,Q,H,Pinf,Pstar, ...
Y,start,Z,kalman_tol,riccati_tol,...
data_index,number_of_observations,...
no_more_missing_observations);
end
end
if isnan(LIK)
cost_flag = 0;
return
end
if imag(LIK)~=0
likelihood = bayestopt_.penalty;
else
likelihood = LIK;
end
% ------------------------------------------------------------------------------
% Adds prior if necessary
% ------------------------------------------------------------------------------
lnprior = priordens(xparam1,bayestopt_.pshape,bayestopt_.p6,bayestopt_.p7,bayestopt_.p3,bayestopt_.p4);
fval = (likelihood-lnprior);
likelihood
lnprior
fval
pause
LIKDLL=logposterior(xparam1,Y,mexext)
pause
options_.kalman_algo = kalman_algo;

View File

@ -0,0 +1,967 @@
C =[
-7.4073
-6.1860
-6.5983
-5.6088
-5.0547
-4.4774
-3.8081
-3.8425
-2.4178
-1.9835
-1.0395
-0.1583
-0.0397
0.3505
-0.1879
-0.0067
0.0478
-1.2247
-1.4349
-0.7973
-0.0461
0.5844
1.1372
1.3801
1.8023
2.2972
2.0469
2.5435
2.8169
3.2007
2.6705
3.0518
3.2445
3.8443
3.8525
4.9494
4.2770
4.9532
5.1441
3.7124
3.9880
3.6926
2.6005
1.8679
1.9085
1.5563
1.2308
0.3264
-0.2208
-0.2483
-0.4082
-1.0315
-1.6030
-1.5499
-1.3777
-2.1675
-2.5138
-2.8820
-2.6958
-2.4719
-1.9854
-1.7954
-2.2362
-1.0595
-0.8808
-0.8548
-1.2839
-0.1363
0.2104
0.8810
0.3555
0.4766
1.3269
1.4506
1.4308
1.6263
1.9842
2.3948
2.8710
3.0177
2.9305
3.1739
3.7380
3.8285
3.3342
3.7447
3.7830
3.1039
2.8413
3.0338
0.3669
0.0847
0.0104
0.2115
-0.6649
-0.9625
-0.7330
-0.8664
-1.4441
-1.0179
-1.2729
-1.9539
-1.4427
-2.0371
-1.9764
-2.5654
-2.8570
-2.5842
-3.0427
-2.8312
-2.3320
-2.2768
-2.1816
-2.1043
-1.8969
-2.2388
-2.1679
-2.1172
];
E =[
0.6263
0.7368
0.7477
1.0150
0.6934
0.4135
0.3845
0.2380
0.2853
0.5999
0.8622
1.2116
1.4921
1.5816
1.7259
1.6276
1.2422
0.8084
0.4710
-0.3704
-0.6427
-0.5323
-0.5562
-0.3651
-0.4356
-0.7164
-0.5816
-0.4635
-0.8456
-0.9708
-0.7138
-0.7499
-0.6941
-0.6656
-0.2912
-0.1650
0.0774
0.2307
0.4484
0.4942
0.4653
0.2196
0.1736
-0.1595
-0.3918
-0.4611
-0.8493
-0.7384
-1.0604
-1.2166
-1.7187
-1.6932
-1.7830
-1.7035
-2.2079
-2.3769
-2.2511
-2.1093
-2.4638
-2.4027
-2.1313
-1.9199
-1.7941
-1.4661
-1.2269
-1.0392
-1.0725
-0.7156
-0.4778
-0.4233
-0.0409
0.1620
0.4280
0.5873
1.0323
1.3420
1.6902
2.0680
2.8219
3.2511
3.2930
3.5633
3.8992
3.6874
3.2849
3.1614
2.6221
2.5067
1.9223
1.1777
0.4483
-0.0661
-0.4424
-0.9000
-1.1478
-1.2047
-1.1412
-1.2383
-1.1048
-0.9716
-0.9287
-1.0057
-1.0827
-1.0200
-1.0072
-1.1740
-1.2809
-1.1086
-0.9866
-0.8947
-0.5875
-0.2329
0.1493
0.4906
0.8400
1.0720
1.2648
1.5431
];
I =[
2.6617
2.4325
1.9592
3.2530
2.9949
3.7918
4.7444
4.8289
5.5983
7.8923
9.4297
9.5010
10.0150
10.0413
9.6046
6.4766
5.9647
3.0114
0.5683
-2.1226
-2.1855
-0.8329
-1.5207
-1.3419
-1.7897
-0.1476
0.4675
-1.6516
-1.5419
-1.3050
-1.2451
-0.7815
-0.7796
-0.3612
-2.4072
1.1162
1.1383
3.4132
5.0356
2.8016
2.1734
0.9366
-0.7050
-1.5021
-2.9868
-6.0237
-6.2589
-6.9138
-8.2340
-9.2589
-9.2465
-9.6988
-9.7782
-10.5645
-10.7544
-13.1583
-12.2718
-12.0131
-13.5983
-12.3579
-10.9146
-11.1572
-12.4935
-9.4393
-8.5535
-7.3723
-10.0169
-6.6088
-5.2045
-4.1024
-2.8472
-1.3139
0.0477
1.5629
3.6947
4.0327
4.1320
7.1400
9.1036
8.5609
7.6576
8.8022
8.9611
10.0871
9.4797
9.3964
10.0363
8.6340
6.6522
4.4471
0.2854
-2.1879
-2.9879
-4.1021
-2.7713
-2.2281
-1.2908
-0.3250
0.6534
0.3942
0.3534
-0.1532
-1.7936
0.4909
0.3634
0.4290
-0.9709
0.1942
0.6103
1.4426
2.7225
1.7525
3.2780
3.5985
4.9011
5.3312
6.4402
6.6529
];
L =[
0.6263
0.7368
0.7477
1.0150
0.6934
0.4135
0.3845
0.2380
0.2853
0.5999
0.8622
1.2116
1.4921
1.5816
1.7259
1.6276
1.2422
0.8084
0.4710
-0.3704
-0.6427
-0.5323
-0.5562
-0.3651
-0.4356
-0.7164
-0.5816
-0.4635
-0.8456
-0.9708
-0.7138
-0.7499
-0.6941
-0.6656
-0.2912
-0.1650
0.0774
0.2307
0.4484
0.4942
0.4653
0.2196
0.1736
-0.1595
-0.3918
-0.4611
-0.8493
-0.7384
-1.0604
-1.2166
-1.7187
-1.6932
-1.7830
-1.7035
-2.2079
-2.3769
-2.2511
-2.1093
-2.4638
-2.4027
-2.1313
-1.9199
-1.7941
-1.4661
-1.2269
-1.0392
-1.0725
-0.7156
-0.4778
-0.4233
-0.0409
0.1620
0.4280
0.5873
1.0323
1.3420
1.6902
2.0680
2.8219
3.2511
3.2930
3.5633
3.8992
3.6874
3.2849
3.1614
2.6221
2.5067
1.9223
1.1777
0.4483
-0.0661
-0.4424
-0.9000
-1.1478
-1.2047
-1.1412
-1.2383
-1.1048
-0.9716
-0.9287
-1.0057
-1.0827
-1.0200
-1.0072
-1.1740
-1.2809
-1.1086
-0.9866
-0.8947
-0.5875
-0.2329
0.1493
0.4906
0.8400
1.0720
1.2648
1.5431
];
PIE =[
-1.0113
-0.8305
0.2332
-0.8746
-0.7978
-0.9220
-0.2487
-0.7749
-0.5460
-0.5347
0.5050
-0.0334
0.6756
0.8791
0.7267
1.0997
1.1750
1.1927
0.4420
0.5357
0.0345
0.0196
0.3371
0.9379
1.2160
0.3393
0.5813
0.7410
0.3374
0.2616
0.4025
0.4799
0.5981
-0.1523
0.4458
0.2182
0.9793
0.7562
1.0064
0.8203
0.6966
0.3352
0.6581
0.6111
0.9833
1.1991
0.9562
0.3868
0.2939
0.2471
0.8331
0.0715
0.3910
0.3301
0.2547
-0.2702
-0.2998
-0.1953
-0.2293
-0.3284
0.0480
-0.0374
0.3253
-0.3434
-0.3892
-0.7178
-0.4758
-0.6794
-0.8505
-0.3512
-0.4436
-0.5101
-0.4574
-0.2696
-0.1047
-0.5745
-0.2989
-0.0063
0.0088
-0.1184
-0.1506
-0.4073
0.2674
0.2896
0.0669
0.1166
-0.1699
-0.2518
-0.0562
-0.3269
-0.0703
-0.1046
-0.4888
-0.3524
-0.2485
-0.5870
-0.4546
-0.3970
-0.2353
-0.0352
-0.2171
-0.3754
-0.4322
-0.4572
-0.4903
-0.4518
-0.6435
-0.6304
-0.4148
-0.2892
-0.4318
-0.6010
-0.4148
-0.4315
-0.3531
-0.8053
-0.4680
-0.4263
];
R =[
-1.0750
-1.1540
-1.3682
-1.4569
-1.3490
-1.4011
-1.6486
-1.6968
-1.6976
-1.2567
-1.1392
-0.7783
-0.3021
-0.0435
0.0066
-0.0043
0.1029
-0.0628
-0.5358
-0.9627
-1.1079
-1.0918
-0.9966
-0.6223
-0.3616
-0.2711
-0.0997
-0.2810
-0.3710
-0.3167
-0.5301
-0.5826
-0.3194
-0.2713
-0.5287
-0.2432
0.1098
0.5349
0.7094
0.8415
0.6226
0.7376
0.9316
1.4370
1.5853
1.4267
1.1783
1.2046
0.9689
0.7918
0.6315
0.5950
0.6853
0.7171
0.5887
0.4873
0.4027
0.3489
0.2934
0.3060
0.1741
0.0348
0.0771
-0.1005
-0.1518
-0.1104
-0.0681
-0.0059
0.0256
0.0404
-0.1721
-0.2002
0.0015
0.1249
0.3738
0.4320
0.5579
0.8186
0.8727
0.7356
0.7243
0.8635
0.9058
0.7656
0.7936
0.8631
0.9074
0.9547
1.2045
1.0850
0.9178
0.5242
0.3178
0.1472
0.0227
-0.0799
-0.0611
-0.0140
0.1132
0.1774
0.0782
0.0436
-0.1596
-0.2691
-0.2895
-0.3791
-0.4020
-0.4166
-0.4037
-0.3636
-0.4075
-0.4311
-0.4470
-0.5111
-0.6274
-0.7261
-0.6974
-0.5012
];
W =[
-14.8791
-13.2300
-13.5037
-13.0249
-11.2546
-10.0148
-8.8586
-8.5739
-7.7851
-6.7136
-5.5878
-4.6881
-3.8039
-3.0366
-2.7342
-1.3135
-0.7387
-0.1131
-0.2769
0.8696
1.8855
2.3667
2.4942
3.2049
3.9682
5.1500
4.7047
4.7827
5.3377
5.6614
5.2813
5.2967
5.5175
6.1526
5.6627
6.0694
6.5824
6.9032
6.7849
6.6896
6.6201
6.9933
5.8959
6.7419
6.9999
6.4009
5.5083
5.1054
5.2813
4.5790
3.9589
3.8599
3.8978
2.7957
3.2480
1.4634
1.9219
1.8398
1.9279
1.8316
1.6092
1.2741
0.2031
-0.0236
-0.1004
-0.3034
-1.0273
-0.2205
0.0458
0.2386
-0.0977
-0.3145
-0.1416
-0.7009
-0.9082
-0.8802
-0.5644
-0.5852
-0.5346
0.0652
0.1301
0.3444
-0.3592
0.8096
0.9644
1.0289
1.2781
1.2298
2.2134
2.0808
0.4925
0.6506
0.5531
0.2456
-0.5351
-0.8183
-0.8967
-0.7268
-1.0738
-1.2844
-1.4338
-1.6995
-1.7085
-2.2889
-2.1018
-2.4273
-2.4609
-2.1407
-2.3847
-3.1689
-4.5581
-4.1027
-4.2436
-4.8836
-5.9660
-4.9971
-5.2386
-5.6618
];
Y =[
-4.9347
-4.6205
-5.2198
-4.5937
-3.8015
-3.6643
-2.7239
-2.7524
-2.0634
-1.0112
0.0530
0.7623
1.7927
2.1486
2.4866
2.1456
2.1671
-0.0254
-1.6716
-1.9673
-1.6109
-1.0292
-0.1222
0.7329
1.1234
2.0603
1.7998
1.4820
1.1732
1.6424
1.5382
2.1399
2.0127
2.7210
2.4966
3.5249
3.6237
4.2011
4.5634
3.3442
2.7761
1.9812
1.3779
1.4616
1.3029
0.7594
0.3695
0.0832
-0.8118
-1.4557
-1.4850
-1.2346
-1.5696
-1.3785
-0.7682
-2.0308
-1.7778
-1.7801
-2.1711
-1.7469
-1.3413
-1.3352
-2.4390
-1.2125
-1.1695
-1.0891
-2.4753
-1.3503
-0.9412
-0.1470
0.0026
0.1108
0.6890
1.3520
1.6018
2.0667
1.7625
2.6658
3.4048
3.2507
3.4251
3.2174
3.1903
3.3396
3.1358
2.8625
3.3546
2.4609
1.9534
0.9962
-0.7904
-1.1672
-1.2586
-1.3593
-1.3443
-0.9413
-0.6023
-0.4516
-0.5129
-0.8741
-1.0784
-1.4091
-1.3627
-1.5731
-1.6037
-1.8814
-2.1482
-1.3597
-1.1855
-1.1122
-0.8424
-0.9747
-1.1385
-1.4548
-1.4284
-1.4633
-1.0621
-0.7871
];

View File

@ -0,0 +1,182 @@
//options_.usePartInfo=1;
var MC E EF R_KF QF CF IF YF LF PIEF WF RF R_K Q C I Y L PIE W R EE_A PIE_BAR EE_B EE_G EE_L EE_I KF K one BIGTHETA;
varexo E_A E_B E_G E_L E_I ETA_R E_PIE_BAR ETA_Q ETA_P ETA_W ;
parameters
xi_e lambda_w alpha czcap beta phi_i tau sig_c hab ccs cinvs phi_y gamma_w xi_w gamma_p xi_p sig_l r_dpi
r_pie r_dy r_y rho rho_a rho_pb rho_b rho_g rho_l rho_i LMP ;
alpha=.30;
beta=.99;
tau=0.025;
ccs=0.6;
cinvs=.22; //% alpha*(tau+ctrend)/R_K R_K=ctrend/beta-1+tau
lambda_w = 0.5;
phi_i= 6.771;
sig_c= 1.353;
hab= 0.573;
xi_w= 0.737;
sig_l= 2.400;
xi_p= 0.908;
xi_e= 0.599;
gamma_w= 0.763;
gamma_p= 0.469;
czcap= 0.169;
phi_y= 1.408;
r_pie= 1.684;
r_dpi= 0.14;
rho= 0.961;
r_y= 0.099;
r_dy= 0.159;
rho_a= 0.823;
rho_b= 0.855;
rho_g= 0.949;
rho_l= 0.889;
rho_i= 0.927;
rho_pb= 0.924;
LMP = 0.0 ; //NEW.
model(linear, use_dll);
CF = (1/(1+hab))*(CF(1)+hab*CF(-1))-((1-hab)/((1+hab)*sig_c))*(RF-PIEF(1)-EE_B) ;
0 = alpha*R_KF+(1-alpha)*WF -EE_A ;
PIEF = 0*one;
IF = (1/(1+beta))* (( IF(-1) + beta*(IF(1)))+(1/phi_i)*QF)+0*ETA_Q+EE_I ;
QF = -(RF-PIEF(1))+(1-beta*(1-tau))*((0+czcap)/czcap)*R_KF(1)+beta*(1-tau)*QF(1) +0*EE_I ;
KF = (1-tau)*KF(-1)+tau*IF(-1) ;
YF = (ccs*CF+cinvs*IF)+EE_G ;
YF = 1*phi_y*( alpha*KF+alpha*(1/czcap)*R_KF+(1-alpha)*LF+EE_A ) ;
WF = (sig_c/(1-hab))*(CF-hab*CF(-1)) + sig_l*LF - EE_L ;
LF = R_KF*((1+czcap)/czcap)-WF+KF ;
EF = EF(-1)+EF(1)-EF+(LF-EF)*((1-xi_e)*(1-xi_e*beta)/(xi_e));
C = (hab/(1+hab))*C(-1)+(1/(1+hab))*C(1)-((1-hab)/((1+hab)*sig_c))*(R-PIE(1)-EE_B) ;
I = (1/(1+beta))* (( I(-1) + beta*(I(1)))+(1/phi_i)*Q )+1*ETA_Q+1*EE_I ;
Q = -(R-PIE(1))+(1-beta*(1-tau))*((0+czcap)/czcap)*R_K(1)+beta*(1-tau)*Q(1) +EE_I*0+0*ETA_Q ;
K = (1-tau)*K(-1)+tau*I(-1) ;
Y = (ccs*C+cinvs*I)+ EE_G ;
Y = phi_y*( alpha*K+alpha*(1/czcap)*R_K+(1-alpha)*L ) +phi_y*EE_A ;
PIE = (1/(1+beta*gamma_p))*
(
(beta)*(PIE(1)) +(gamma_p)*(PIE(-1))
+((1-xi_p)*(1-beta*xi_p)/(xi_p))*(MC)
) + ETA_P ;
MC = alpha*R_K+(1-alpha)*W -EE_A;
W = (1/(1+beta))*(beta*W(+1)+W(-1))
+(beta/(1+beta))*(PIE(+1))
-((1+beta*gamma_w)/(1+beta))*(PIE)
+(gamma_w/(1+beta))*(PIE(-1))
-(1/(1+beta))*(((1-beta*xi_w)*(1-xi_w))/(((1+(((1+lambda_w)*sig_l)/(lambda_w))))*xi_w))*(W-sig_l*L-(sig_c/(1-hab))*(C-hab*C(-1))+EE_L)
+ETA_W;
L = R_K*((1+czcap)/czcap)-W+K ;
// R = r_dpi*(PIE-PIE(-1))
// +(1-rho)*(r_pie*(PIE(-1)-PIE_BAR)+r_y*(Y-YF))
// +r_dy*(Y-YF-(Y(-1)-YF(-1)))
// +rho*(R(-1)-PIE_BAR)
// +PIE_BAR
// +ETA_R;
R =
r_dpi*(PIE-PIE(-1))
+(1-rho)*(r_pie*(BIGTHETA)+r_y*(Y-YF))
+r_dy*(Y-YF-(Y(-1)-YF(-1)))
+rho*(R(-1)-PIE_BAR)
+PIE_BAR
+ETA_R;
E = E(-1)+E(1)-E+(L-E)*((1-xi_e)*(1-xi_e*beta)/(xi_e));
EE_A = (rho_a)*EE_A(-1) + E_A;
PIE_BAR = rho_pb*PIE_BAR(-1)+ E_PIE_BAR ;
EE_B = rho_b*EE_B(-1) + E_B ;
EE_G = rho_g*EE_G(-1) + E_G ;
EE_L = rho_l*EE_L(-1) + E_L ;
EE_I = rho_i*EE_I(-1) + E_I ;
one = 0*one(-1) ;
LMP*BIGTHETA(1) = BIGTHETA - (PIE(-1) - PIE_BAR) ;
end;
shocks;
var E_A; stderr 0.598;
var E_B; stderr 0.336;
var E_G; stderr 0.325;
var E_I; stderr 0.085;
var E_L; stderr 3.520;
var ETA_P; stderr 0.160;
var ETA_W; stderr 0.289;
var ETA_R; stderr 0.081;
var ETA_Q; stderr 0.604;
var E_PIE_BAR; stderr 0.017;
end;
//stoch_simul(irf=20) Y C PIE R W R_K L Q I K ;
// stoch_simul generates what kind of standard errors for the shocks ?
//steady;
//check;
//stoch_simul(periods=200,irf=20,simul_seed=3) Y C PIE MC R W R_K E L I ;
//datatomfile('ddd',[]);
// new syntax
estimated_params;
// PARAM NAME, INITVAL, LB, UB, PRIOR_SHAPE, PRIOR_P1, PRIOR_P2, PRIOR_P3, PRIOR_P4, JSCALE
// PRIOR_SHAPE: BETA_PDF, GAMMA_PDF, NORMAL_PDF, INV_GAMMA_PDF
stderr E_A,0.543,0.01,4,INV_GAMMA_PDF,0.4,2;
stderr E_PIE_BAR,0.072,0.001,4,INV_GAMMA_PDF,0.02,10;
stderr E_B,0.2694,0.01,4,INV_GAMMA_PDF,0.2,2;
stderr E_G,0.3052,0.01,4,INV_GAMMA_PDF,0.3,2;
stderr E_L,1.4575,0.1,6,INV_GAMMA_PDF,1,2;
stderr E_I,0.1318,0.01,4,INV_GAMMA_PDF,0.1,2;
stderr ETA_R,0.1363,0.01,4,INV_GAMMA_PDF,0.1,2;
stderr ETA_Q,0.4842,0.01,4,INV_GAMMA_PDF,0.4,2;
stderr ETA_P,0.1731,0.01,4,INV_GAMMA_PDF,0.15,2;
stderr ETA_W,0.2462,0.1,4,INV_GAMMA_PDF,0.25,2;
rho_a,.9722,.1,.9999,BETA_PDF,0.85,0.1;
rho_pb,.85,.1,.999,BETA_PDF,0.85,0.1;
rho_b,.7647,.1,.99,BETA_PDF,0.85,0.1;
rho_g,.9502,.1,.9999,BETA_PDF,0.85,0.1;
rho_l,.9542,.1,.9999,BETA_PDF,0.85,0.1;
rho_i,.6705,.1,.99,BETA_PDF,0.85,0.1;
phi_i,5.2083,1,15,NORMAL_PDF,4,1.5;
sig_c,0.9817,0.25,3,NORMAL_PDF,1,0.375;
hab,0.5612,0.3,0.95,BETA_PDF,0.7,0.1;
xi_w,0.7661,0.3,0.9,BETA_PDF,0.75,0.05;
sig_l,1.7526,0.5,5,NORMAL_PDF,2,0.75;
xi_p,0.8684,0.3,0.95,BETA_PDF,0.75,0.05;
xi_e,0.5724,0.1,0.95,BETA_PDF,0.5,0.15;
gamma_w,0.6202,0.1,0.99,BETA_PDF,0.75,0.15;
gamma_p,0.6638,0.1,0.99,BETA_PDF,0.75,0.15;
czcap,0.2516,0.01,2,NORMAL_PDF,0.2,0.075;
phi_y,1.3011,1.001,2,NORMAL_PDF,1.45,0.125;
r_pie,1.4616,1.2,2,NORMAL_PDF,1.7,0.1;
r_dpi,0.1144,0.01,0.5,NORMAL_PDF,0.3,0.1;
rho,0.8865,0.5,0.99,BETA_PDF,0.8,0.10;
r_y,0.0571,0.01,0.2,NORMAL_PDF,0.125,0.05;
r_dy,0.2228,0.05,0.5,NORMAL_PDF,0.0625,0.05;
end;
varobs Y C I E PIE W R;
//estimation(datafile=rawdata_euromodel_1,presample=40, first_obs=1, nobs=118, lik_init=2, mode_compute=1,mh_replic=0);
estimation(datafile=rawdata_euromodel_1,presample=40, first_obs=1, nobs=118,mh_jscale=0.2,mh_replic=150000, mode_check); //
//stoch_simul(periods=200,irf=20,simul_seed=3) Y C PIE R W R_K L Q I K ;