Commit 05b9c70a authored by Michel Juillard's avatar Michel Juillard
Browse files

removing command CALIB: it has not been working for a long time. Calibration...

removing command CALIB: it has not been working for a long time. Calibration is in fact a special case of method of moments and should be replaced by such a method.
parent 256ff761
function M_.Sigma_e = calib(var_indices,targets,var_weights,nar,cova,M_.Sigma_e)
% Copyright (C) 2005 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 oo_ M_ vx
ncstr = 0;
ni = size(var_indices,1);
for i=1:nar+3
ncstr = ncstr + size(var_indices{i},1);
end
if cova
if ncstr < M_.exo_nbr*(M_.exo_nbr+1)/2
error(['number of preset variances is smaller than number of shock' ...
' variances and covariances to be estimated !'])
end
else
if ncstr < M_.exo_nbr
error(['number of preset variances is smaller than number of shock' ...
' variances to be estimated !'])
end
end
% computes approximate solution at order 1
dr = resol(oo_.steady_state,0,0,1);
ghx = dr.ghx;
ghu = dr.ghu;
npred = dr.npred;
nstatic = dr.nstatic;
kstate = dr.kstate;
order = dr.order_var;
iv(order) = [1:M_.endo_nbr];
iv = iv';
nx = size(ghx,2);
ikx = [nstatic+1:nstatic+npred];
A = zeros(nx,nx);
A(1:npred,:)=ghx(ikx,:);
offset_r = npred;
offset_c = 0;
i0 = find(kstate(:,2) == M_.maximum_lag+1);
n0 = size(i0,1);
for i=M_.maximum_lag:-1:2
i1 = find(kstate(:,2) == i);
n1 = size(i1,1);
j = zeros(n1,1);
for j1 = 1:n1
j(j1) = find(kstate(i0,1)==kstate(i1(j1),1));
end
A(offset_r+1:offset_r+n1,offset_c+j)=eye(n1);
offset_r = offset_r + n1;
offset_c = offset_c + n0;
i0 = i1;
n0 = n1;
end
ghu1 = [ghu(ikx,:);zeros(nx-npred,M_.exo_nbr)];
% IA = speye(nx*nx)-kron(A,A);
% kron_ghu = kron(ghu1,ghu1);
% vx1 such that vec(sigma_x) = vx1 * vec(M_.Sigma_e) (predetermined vars)
vx1 = [];
% vx1 = IA\kron_ghu;
IA = [];
kron_ghu = [];
% computes required variables and indices among required variables
iiy = [];
for i=1:nar+3
if i ~= 3 & ~isempty(var_indices{i})
iiy = union(iiy, iv(var_indices{i}(:,1)));
end
end
if ~isempty(var_indices{2})
iiy = union(iiy, iv(var_indices{2}(:,2)));
end
ny = size(iiy,1);
for i=1:nar+3
if i ~= 3 & ~isempty(var_indices{i})
var_indices{i}(:,1) = indnv(iv(var_indices{i}(:,1)),iiy);
end
if i ~= 2 & i ~= 3 & ~isempty(var_indices{i})
var_indices{i} = sub2ind([ny ny],var_indices{i},var_indices{i});
end
end
if ~isempty(var_indices{2})
var_indices{2}(:,2) = indnv(iv(var_indices{2}(:,2)),iiy);
var_indices{2} = sub2ind([ny ny],var_indices{2}(:,1),var_indices{2}(:,2));
end
if ~isempty(var_indices{3})
var_indices{3} = sub2ind([M_.exo_nbr M_.exo_nbr],var_indices{3}(:,1),var_indices{3}(:,2));
end
if isempty(M_.Sigma_e)
M_.Sigma_e = 0.01*eye(M_.exo_nbr);
b = 0.1*ghu1*ghu1';
else
b = ghu1*M_.Sigma_e*ghu1';
M_.Sigma_e = chol(M_.Sigma_e+1e-14*eye(M_.exo_nbr));
end
options=optimset('LargeScale','on','MaxFunEvals',20000*ny,'TolX',1e-4, ...
'TolFun',1e-4,'Display','Iter','MaxIter',10000);
% [M_.Sigma_e,f]=fminunc(@calib_obj,M_.Sigma_e,options,A,ghu1,ghx(iiy,:),ghu(iiy,:),targets,var_weights,var_indices,nar);
[M_.Sigma_e,f]=fmincon(@calib_obj,diag(M_.Sigma_e).^2,-eye(M_.exo_nbr),zeros(M_.exo_nbr,1),[],[],[],[],[],options,A,ghu1,ghx(iiy,:),ghu(iiy,:),targets,var_weights,var_indices,nar);
M_.Sigma_e = diag(M_.Sigma_e);
objective = calib_obj2(diag(M_.Sigma_e),A,ghu1,ghx(iiy,:),ghu(iiy,:),targets,var_weights,var_indices,nar);
disp('CALIBRATION')
disp('')
if ~isempty(var_indices{1})
disp(sprintf('%16s %14s %14s %14s %14s','Std. Dev','Target','Obtained','Diff'));
str = ' ';
for i=1:size(var_indices{1},1)
[i1,i2] = ind2sub([ny ny],var_indices{1}(i));
str = sprintf('%16s: %14.2f %14.2f %14.2f',M_.endo_names(order(iiy(i1)),:),targets{1}(i),objective{1}(i),objective{1}(i)-targets{1}(i));
disp(str);
end
end
if ~isempty(var_indices{2})
disp(sprintf('%32s %14s %14s','Correlations','Target','Obtained','Diff'));
str = ' ';
for i=1:size(var_indices{2},1)
[i1,i2]=ind2sub([ny ny],var_indices{2}(i));
str = sprintf('%16s,%16s: %14.2f %14.2f %14.2f',M_.endo_names(order(iiy(i1)),:), ...
M_.endo_names(order(iiy(i2)),:),targets{2}(i),objective{2}(i),objective{2}(i)-targets{2}(i));
disp(str);
end
end
if ~isempty(var_indices{3})
disp(sprintf('%32s %16s %16s','Constrained shocks (co)variances','Target','Obtained'));
str = ' ';
for i=1:size(var_indices{3},1)
[i1,i2]=ind2sub([M_.exo_nbr M_.exo_nbr],var_indices{3}(i));
if i1 == i2
str = sprintf('%32s: %16.4f %16.4f',M_.exo_name(order(i1),:), ...
targets{3}(i),objective{3}(i));
else
str = sprintf('%16s,%16s: %16.4f %16.4f',M_.exo_name(order(i1),:), ...
M_.exo_name(order(i2), :),targets{3}(i),objective{3}(i));
end
disp(str);
end
end
flag = 1;
for j=4:nar+3
if ~isempty(var_indices{j})
if flag
disp(sprintf('%16s %16s %16s','Autocorrelations','Target','Obtained'));
str = ' ';
flag = 0;
end
for i=1:size(var_indices{j},1)
[i1,i2] = ind2sub([ny ny],var_indices{j}(i));
str = sprintf('%16s(%d): %16.4f %16.4f',M_.endo_names(order(iiy(i1)),:), ...
j-3,targets{j}(i),objective{j}(i));
disp(str);
end
end
end
disp('');
disp('Calibrated variances')
str = ' ';
for i=1:M_.exo_nbr
str = [str sprintf('%16s',M_.exo_name(i,:))];
end
disp(str);
disp('');
str = ' ';
for i=1:M_.exo_nbr
str = [str sprintf('%16f',M_.Sigma_e(i,i))];
end
disp(str);
% 10/9/02 MJ
\ No newline at end of file
function f=calib_obj(M_.Sigma_e,A,ghu1,ghx,ghu,targets,var_weights,iy,nar)
% targets and iy order: 1) variances 2) correlations
% 3) constraints on M_.Sigma_e itself 4) autocorrelations
% Copyright (C) 2005-2008 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 vx fold options_
oo_.gamma_y = cell(nar+1,1);
% M_.Sigma_e = M_.Sigma_e'*M_.Sigma_e;
M_.Sigma_e=diag(M_.Sigma_e);
nx = size(ghx,2);
b=ghu1*M_.Sigma_e*ghu1';
vx = [];
if isempty(vx)
vx = lyapunov_symm(A,b,options_.qz_criterium,options_.lyapunov_complex_threshold);
else
[vx,status] = bicgstab_(@f_var,b(:),vx(:),1e-8,50,A,nx);
if status
vx = lyapunov_symm(A,b,options_.qz_criterium,options_.lyapunov_complex_threshold);
else
vx=reshape(vx,nx,nx);
end
end
oo_.gamma_y{1} = ghx*vx*ghx'+ ghu*M_.Sigma_e*ghu';
f = 0;
if ~isempty(targets{1})
e = targets{1}-sqrt(oo_.gamma_y{1}(iy{1}));
f = e'*(var_weights{1}.*e);
end
sy = sqrt(diag(oo_.gamma_y{1}));
sy = sy *sy';
if ~isempty(targets{2})
e = targets{2}-oo_.gamma_y{1}(iy{2})./(sy(iy{2})+1e-10);
f = f+e'*(var_weights{2}.*e);
end
if ~isempty(targets{3})
e = targets{3}-sqrt(M_.Sigma_e(iy{3}));
f = f+e'*(var_weights{3}.*e);
end
% autocorrelations
if nar > 0
vxy = (A*vx*ghx'+ghu1*M_.Sigma_e*ghu');
oo_.gamma_y{2} = ghx*vxy./(sy+1e-10);
if ~isempty(targets{4})
e = targets{4}-oo_.gamma_y{2}(iy{4});
f = f+e'*(var_weights{4}.*e);
end
for i=2:nar
vxy = A*vxy;
oo_.gamma_y{i+1} = ghx*vxy./(sy+1e-10);
if ~isempty(targets{i+3})
e = targets{i+3}-oo_.gamma_y{i+1}(iy{i+3});
f = f+e'*(var_weights{i+3}.*e);
end
end
end
if isempty(fold) | f < 2*fold
fold = f;
vxold = vx;
end
% 11/04/02 MJ generalized for correlations, autocorrelations and
% constraints on M_.Sigma_e
% 01/25/03 MJ targets std. dev. instead of variances
function objective=calib_obj2(M_.Sigma_e,A,ghu1,ghx,ghu,targets,var_weights,iy,nar)
% targets and iy order: 1) variances 2) correlations
% 3) constraints on M_.Sigma_e itself 4) autocorrelations
% Copyright (C) 2005-2008 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 vx fold options_
objective = cell (nar+3);
oo_.gamma_y = cell(nar+1,1);
M_.Sigma_e=diag(M_.Sigma_e);
nx = size(ghx,2);
b=ghu1*M_.Sigma_e*ghu1';
vx = lyapunov_symm(A,b,options_.qz_criterium,options_.lyapunov_complex_threshold);
oo_.gamma_y{1} = ghx*vx*ghx'+ ghu*M_.Sigma_e*ghu';
if ~isempty(targets{1})
objective{1} = sqrt(oo_.gamma_y{1}(iy{1}));
end
sy = sqrt(diag(oo_.gamma_y{1}));
sy = sy *sy';
if ~isempty(targets{2})
objective{2} = oo_.gamma_y{1}(iy{2})./(sy(iy{2})+1e-10);
end
if ~isempty(targets{3})
objective{3} = M_.Sigma_e(iy{3});
end
% autocorrelations
if nar > 0
vxy = (A*vx*ghx'+ghu1*M_.Sigma_e*ghu');
oo_.gamma_y{2} = ghx*vxy./(sy+1e-10);
if ~isempty(targets{4})
objective{4} = oo_.gamma_y{2}(iy{4});
end
for i=2:nar
vxy = A*vxy;
oo_.gamma_y{i+1} = ghx*vxy./(sy+1e-10);
if ~isempty(targets{i+3})
objecitve{i+3} = oo_.gamma_y{i+1}(iy{i+3});
end
end
end
% 11/04/02 MJ generalized for correlations, autocorrelations and
% constraints on M_.Sigma_e
\ No newline at end of file
......@@ -737,7 +737,7 @@ CalibStatement::CalibStatement(int covar_arg) : covar(covar_arg)
void
CalibStatement::writeOutput(ostream &output, const string &basename) const
{
output << "M_.Sigma_e=calib(calib_var_index,calib_targets,calib_weights," << covar << ",Sigma_e_);\n";
output << "M_.Sigma_e=calib(calib_var_index,calib_targets,calib_weights,options_.ar," << covar << ",M_.Sigma_e);\n";
}
OsrParamsStatement::OsrParamsStatement(const SymbolList &symbol_list_arg) :
......
......@@ -93,7 +93,7 @@ class ParsingDriver;
%token BVAR_PRIOR_DECAY BVAR_PRIOR_FLAT BVAR_PRIOR_LAMBDA
%token BVAR_PRIOR_MU BVAR_PRIOR_OMEGA BVAR_PRIOR_TAU BVAR_PRIOR_TRAIN
%token BVAR_REPLIC BYTECODE
%token CALIB CALIB_VAR CHANGE_TYPE CHECK CONDITIONAL_FORECAST CONDITIONAL_FORECAST_PATHS CONF_SIG CONSTANT CONTROLLED_VAREXO CORR COVAR CUTOFF
%token CHANGE_TYPE CHECK CONDITIONAL_FORECAST CONDITIONAL_FORECAST_PATHS CONF_SIG CONSTANT CONTROLLED_VAREXO CORR COVAR CUTOFF
%token DATAFILE DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE
%token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT
%token FILENAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS
......@@ -166,7 +166,7 @@ class ParsingDriver;
%type <string_val> non_negative_number signed_number signed_integer
%type <string_val> filename symbol expectation_input
%type <string_val> vec_value_1 vec_value
%type <string_val> calib_arg2 range prior
%type <string_val> range prior
%type <symbol_type_val> change_type_arg
%type <vector_string_val> change_type_var_list
%type <vector_int_val> vec_int_elem vec_int_1 vec_int vec_int_number
......@@ -211,8 +211,6 @@ statement : parameters
| optim_weights
| osr_params
| osr
| calib_var
| calib
| dynatype
| dynasave
| model_comparison
......@@ -1230,31 +1228,6 @@ osr : OSR ';'
{driver.run_osr(); }
;
calib_var : CALIB_VAR ';' calib_var_list END ';' { driver.run_calib_var(); };
calib_var_list : calib_var_list calib_arg1
| calib_arg1
;
calib_arg1 : symbol calib_arg2 EQUAL expression ';'
{ driver.set_calib_var($1, $2, $4); }
| symbol COMMA symbol calib_arg2 EQUAL expression ';'
{ driver.set_calib_covar($1, $3, $4, $6); }
| AUTOCORR symbol '(' INT_NUMBER ')' calib_arg2 EQUAL expression ';'
{ driver.set_calib_ac($2, $4, $6, $8); }
;
calib_arg2 : { $$ = new string("1"); }
| '(' non_negative_number ')'
{ $$ = $2; }
;
calib : CALIB ';'
{ driver.run_calib(0); }
| CALIB '(' COVAR ')' ';'
{ driver.run_calib(1); }
;
dynatype : DYNATYPE '(' filename ')' ';'
{ driver.run_dynatype($3); }
| DYNATYPE '(' filename ')' symbol_list ';'
......
......@@ -128,7 +128,6 @@ string eofbuff;
<INITIAL>stoch_simul {BEGIN DYNARE_STATEMENT; return token::STOCH_SIMUL;}
<INITIAL>dsample {BEGIN DYNARE_STATEMENT; return token::DSAMPLE;}
<INITIAL>Sigma_e {BEGIN DYNARE_STATEMENT; sigma_e = 1; return token::SIGMA_E;}
<INITIAL>calib {BEGIN DYNARE_STATEMENT; return token::CALIB;}
<INITIAL>planner_objective {BEGIN DYNARE_STATEMENT; return token::PLANNER_OBJECTIVE;}
<INITIAL>ramsey_policy {BEGIN DYNARE_STATEMENT; return token::RAMSEY_POLICY;}
<INITIAL>identification {BEGIN DYNARE_STATEMENT; return token::IDENTIFICATION;}
......@@ -171,7 +170,6 @@ string eofbuff;
<INITIAL>estimated_params_bounds {BEGIN DYNARE_BLOCK; return token::ESTIMATED_PARAMS_BOUNDS;}
<INITIAL>observation_trends {BEGIN DYNARE_BLOCK; return token::OBSERVATION_TRENDS;}
<INITIAL>optim_weights {BEGIN DYNARE_BLOCK; return token::OPTIM_WEIGHTS;}
<INITIAL>calib_var {BEGIN DYNARE_BLOCK; return token::CALIB_VAR;}
<INITIAL>homotopy_setup {BEGIN DYNARE_BLOCK; return token::HOMOTOPY_SETUP;}
<INITIAL>conditional_forecast_paths {BEGIN DYNARE_BLOCK; return token::CONDITIONAL_FORECAST_PATHS;}
<INITIAL>svar_identification {BEGIN DYNARE_BLOCK; return token::SVAR_IDENTIFICATION;}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment