v4: changes around restricted state space representation

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@918 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
michel 2006-09-17 15:23:45 +00:00
parent 8fcac85ea1
commit cb9fb1b353
5 changed files with 82 additions and 85 deletions

View File

@ -90,8 +90,9 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
%------------------------------------------------------------------------------
% 2. call model setup & reduction program
%------------------------------------------------------------------------------
[T,R,SteadyState,info] = dynare_resolve;
rs = bayestopt_.restrict_state;
[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;
@ -101,8 +102,6 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
cost_flag = 0;
return
end
T = T(rs,rs);
R = R(rs,:);
bayestopt_.mf = bayestopt_.mf1;
if ~options_.noconstant
if options_.loglinear == 1
@ -129,21 +128,19 @@ function [fval,cost_flag,ys,trend_coeff,info] = DsgeLikelihood(xparam1,gend,data
% 3. Initial condition of the Kalman filter
%------------------------------------------------------------------------------
if options_.lik_init == 1 % Kalman filter
Pstar = lyapunov_symm(T,R*Q*transpose(R));
Pstar = lyapunov_symm(T,R*Q*R');
Pinf = [];
elseif options_.lik_init == 2 % Old Diffuse Kalman filter
Pstar = 10*eye(np);
Pinf = [];
elseif options_.lik_init == 3 % Diffuse Kalman filter
Pstar = zeros(np,np);
% ivs = bayestopt_.i_T_var_stable;
% Pstar(ivs,ivs) = lyapunov_symm(T(ivs,ivs),R(ivs,:)*Q* ...
% transpose(R(ivs,:)));
Pstar(ivs,ivs) = lyapunov_symm(T(ivs,ivs),R(ivs,:)*Q* ...
transpose(R(ivs,:)));
Pinf = bayestopt_.Pinf;
ivs = bayestopt_.restrict_var_list_stationary;
R1 = R(ivs,:);
Pstar(ivs,ivs) = lyapunov_symm(T(ivs,ivs),R1*Q*R1');
% Pinf = bayestopt_.Pinf;
% by M. Ratto
% RR=T(:,find(~ismember([1:np],ivs)));
RR=T(:,bayestopt_.restrict_var_list_nonstationary);
i=find(abs(RR)>1.e-10);
R0=zeros(size(RR));
R0(i)=sign(RR(i));

View File

@ -142,63 +142,15 @@ np = estim_params_.np ;
nx = nvx+nvn+ncx+ncn+np;
dr = set_state_space([]);
%% Initialization with unit-root variables
if ~isempty(options_.unit_root_vars)
n_ur = size(options_.unit_root_vars,1);
i_ur = zeros(n_ur,1);
for i=1:n_ur
i1 = strmatch(deblank(options_.unit_root_vars(i,:)),M_.endo_names(dr.order_var,:),'exact');
if isempty(i1)
error('Undeclared variable in unit_root_vars statement')
end
i_ur(i) = i1;
end
if M_.maximum_lag > 1
l1 = flipud([cumsum(M_.lead_lag_incidence(1:M_.maximum_lag-1,dr.order_var),1);ones(1,M_.endo_nbr)]);
n1 = nnz(l1);
bayestopt_.Pinf = zeros(n1,n1);
l2 = find(l1');
l3 = zeros(M_.endo_nbr,M_.maximum_lag);
l3(i_ur,:) = l1(:,i_ur)';
l3 = l3(:);
i_ur1 = find(l3(l2));
i_stable = ones(M_.endo_nbr,1);
i_stable(i_ur) = zeros(n_ur,1);
i_stable = find(i_stable);
bayestopt_.Pinf(i_ur1,i_ur1) = diag(ones(1,length(i_ur1)));
bayestopt_.i_var_stable = i_stable;
l3 = zeros(M_.endo_nbr,M_.maximum_lag);
l3(i_stable,:) = l1(:,i_stable)';
l3 = l3(:);
bayestopt_.i_T_var_stable = find(l3(l2));
else
n1 = M_.endo_nbr;
bayestopt_.Pinf = zeros(n1,n1);
bayestopt_.Pinf(i_ur,i_ur) = diag(ones(1,length(i_ur)));
l1 = ones(M_.endo_nbr,1);
l1(i_ur,:) = zeros(length(i_ur),1);
bayestopt_.i_T_var_stable = find(l1);
end
options_.lik_init = 3;
end % if ~isempty(options_.unit_root_vars)
if isempty(options_.datafile)
error('ESTIMATION: datafile option is missing')
end
nstatic = dr.nstatic;
npred = dr.npred;
nspred = dr.nspred;
if isempty(options_.varobs)
error('ESTIMATION: VAROBS is missing')
end
%% If jscale isn't specified for an estimated parameter, use
%% global option options_.jscale, set to 0.2, by default
k = find(isnan(bayestopt_.jscale));
bayestopt_.jscale(k) = options_.mh_jscale;
%% Read and demean data
rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
%% Setting resticted state space (observed + predetermined variables)
k = [];
k1 = [];
@ -213,7 +165,7 @@ k2 = [k2;[M_.endo_nbr+(1:dr.nspred-dr.npred)]'];
% set restrict_state to postion of observed + state variables
% in expanded state vector
bayestopt_.restrict_state = k2;
bayestopt_.restrict_var_list = k2;
% set mf1 to positions of observed variables in restricted state vector
% for likelihood computation
[junk,bayestopt_.mf1] = ismember(k,k2);
@ -221,6 +173,55 @@ bayestopt_.restrict_state = k2;
% for filtering and smoothing
bayestopt_.mf2 = k;
bayestopt_.mfys = k1;
[junk,ic] = intersect(k2,nstatic+(1:npred)');
bayestopt_.restrict_columns = [ic length(k2)+(1:nspred-npred)]';
aux = dr.transition_auxiliary_variables;
k = find(aux(:,2) > npred);
aux(:,2) = aux(:,2)+sum(k2 <= nstatic);
bayestopt_.restrict_aux = aux;
%% Initialization with unit-root variables
if ~isempty(options_.unit_root_vars)
n_ur = size(options_.unit_root_vars,1);
i_ur = zeros(n_ur,1);
for i=1:n_ur
i1 = strmatch(deblank(options_.unit_root_vars(i,:)),M_.endo_names(dr.order_var,:),'exact');
if isempty(i1)
error('Undeclared variable in unit_root_vars statement')
end
i_ur(i) = i1;
end
[junk,bayestopt_.restrict_var_list_stationary] = ...
setdiff(bayestopt_.restrict_var_list,i_ur);
[junk,bayestopt_.restrict_var_list_nonstationary] = ...
setdiff(bayestopt_.restrict_var_list,i_ur);
if M_.maximum_lag > 1
l1 = flipud([cumsum(M_.lead_lag_incidence(1:M_.maximum_lag-1,dr.order_var),1);ones(1,M_.endo_nbr)]);
l2 = l1(:,restrict_var_list);
il2 = find(l2' > 0);
l2(il2) = (1:length(il2))';
bayestopt_.restict_var_list_stationary = ...
nonzeros(l2(:,restrict_var_list_stationary));
bayestopt_.restict_var_list_nonstationary = ...
nonzeros(l2(:,restrict_var_list_nonstationary));
end
options_.lik_init = 3;
end % if ~isempty(options_.unit_root_vars)
if isempty(options_.datafile)
error('ESTIMATION: datafile option is missing')
end
%% If jscale isn't specified for an estimated parameter, use
%% global option options_.jscale, set to 0.2, by default
k = find(isnan(bayestopt_.jscale));
bayestopt_.jscale(k) = options_.mh_jscale;
%% Read and demean data
rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range);
options_ = set_default_option(options_,'nobs',size(rawdata,1)-options_.first_obs+1);
gend = options_.nobs;

View File

@ -1,5 +1,5 @@
function [A,B,ys,info] = dynare_resolve()
global oo_
function [A,B,ys,info] = dynare_resolve(iv,ic,aux)
global oo_ M_
[oo_.dr,info] = resol(oo_.steady_state,0);
@ -10,5 +10,15 @@ function [A,B,ys,info] = dynare_resolve()
return
end
[A,B] = kalman_transition_matrix(oo_.dr);
if nargin == 0
endo_nbr = M_.endo_nbr;
iv = (1:endo_nbr)';
ic = [ nstatic+(1:npred) endo_nbr+(1:size(dr.ghx,2)-npred) ]';
aux = oo_.dr.transtion_auxiliary_variables(:,2);
k = find(aux > oo_.dr.npred);
aux = aux + nstatic;
aux(k) = aux(k) + oo_.dr.nfrwd;
end
[A,B] = kalman_transition_matrix(oo_.dr,iv,ic,aux);
ys = oo_.dr.ys;

View File

@ -1,4 +1,4 @@
function [loss,vx,info]=osr_obj(x,i_params,weights);
function [loss,vx,info]=osr_obj(x,i_params,i_var,weights);
% objective function for optimal simple rules (OSR)
global M_ oo_ optimal_Q_ it_
% global ys_ Sigma_e_ endo_nbr exo_nbr optimal_Q_ it_ ykmin_ options_
@ -33,16 +33,7 @@ function [loss,vx,info]=osr_obj(x,i_params,weights);
otherwise
end
[A,B] = kalman_transition_matrix(dr);
[vx,ns_var] = lyapunov_symm(A,B*M_.Sigma_e*B');
endo_nbr = M_.endo_nbr;
i_var = (1:endo_nbr)';
i_var(ns_var) = zeros(length(ns_var),1);
i_var = nonzeros(i_var);
vx = vx(i_var,i_var);
weights = weights(dr.order_var,dr.order_var);
weights = sparse(weights(i_var,i_var));
vx = get_variance_of_endogenous_variables(dr,i_var);
loss = weights(:)'*vx(:);

View File

@ -28,7 +28,7 @@ npred = length(pred_var);
nfwrd = length(fwrd_var);
nstatic = length(stat_var);
order_var = [ stat_var; pred_var; both_var; fwrd_var];
inv_order_var(order_var) = (1:M_.endo_nbr)';
inv_order_var(order_var) = (1:M_.endo_nbr);
% building kmask for z state vector in t+1
if M_.maximum_lag > 0
@ -72,7 +72,7 @@ kstate(M_.maximum_lead*M_.endo_nbr+1:end,4) = kiy((M_.maximum_lead+1)*M_.endo_nb
kstate = kstate(i_kmask,:);
dr.order_var = order_var;
dr.inv_order_var = inv_order_var;
dr.inv_order_var = inv_order_var';
dr.nstatic = nstatic;
dr.npred = npred+nboth;
dr.kstate = kstate;
@ -85,9 +85,10 @@ dr.nsfwrd = sum(kstate(:,2) > M_.maximum_lag+1);
% number of predetermined variables in the state vector
dr.nspred = sum(kstate(:,2) <= M_.maximum_lag+1);
% copmutes column position of auxiliary variables for
% compact transition matrix (only state variables)
aux = zeros(0,2);
k0 = kstate(find(kstate(:,2) <= M_.maximum_lag+1),:);;
offset_col = nstatic;
for i=M_.maximum_lag:-1:2
i1 = find(k0(:,2) == i);
n1 = size(i1,1);
@ -95,9 +96,6 @@ for i=M_.maximum_lag:-1:2
for j1 = 1:n1
j(j1) = find(k0(i0,1)==k0(i1(j1),1));
end
if i == M_.maximum_lag-1
offset_col = dr.nstatic+dr.nfwrd;
end
aux = [aux; offset_col+i0(j)];
i0 = i1;
end