- reversed ordering for lower and upper Cholesky
- added tests - fixed various bugs - modified initial values Still unfinishedtime-shift
Michel Juillard
parent
f4b7840739
commit
591e426110
|
|
@ -283,6 +283,15 @@ options_.ms.print_draws = 1;
|
|||
options_.ms.n_draws=1000;
|
||||
options_.ms.thinning_factor=1;
|
||||
options_.ms.proposal_draws = 100000;
|
||||
|
||||
options_.ms.lower_cholesky = 0;
|
||||
options_.ms.upper_cholesky = 0;
|
||||
options_.ms.Qi = [];
|
||||
options_.ms.Ri = [];
|
||||
options_.ms.draws_nbr_burn_in_1 = 30000;
|
||||
options_.ms.draws_nbr_burn_in_2 = 10000;
|
||||
options_.ms.draws_nbr_mean_var_estimate = 200000;
|
||||
options_.ms.draws_nbr_modified_harmonic_mean = 1000000;
|
||||
options_.ms.thinning_factor = 1;
|
||||
options_.ms.dirichlet_scale = 1;
|
||||
% initialize persistent variables in priordens()
|
||||
priordens([],[],[],[],[],[],1);
|
||||
|
|
|
|||
|
|
@ -50,10 +50,11 @@ endif
|
|||
|
||||
ifdef USE_MICHEL_LAPTOP
|
||||
|
||||
#CC = gcc
|
||||
CC = /opt/intel/Compiler/11.0/074/bin/intel64/icc -no-multibyte-chars
|
||||
CC = gcc
|
||||
#CC = /opt/intel/Compiler/11.0/074/bin/intel64/icc -no-multibyte-chars
|
||||
FC = gfortran
|
||||
CFLAGS = -g -static
|
||||
#CFLAGS = -g
|
||||
|
||||
TAO_DIR = $(WORK_DIR)/utilities/TZCcode
|
||||
|
||||
|
|
@ -65,6 +66,7 @@ INTEL_LIBS = -lguide
|
|||
MKL_BASE_DIR = /opt/intel/Compiler/11.0/074/mkl
|
||||
MKL_LIBS_DIR = $(MKL_BASE_DIR)/lib/em64t
|
||||
MKL_LIBS = -lmkl_lapack -lmkl_em64t
|
||||
|
||||
endif
|
||||
|
||||
|
||||
|
|
@ -173,15 +175,16 @@ ERROR_DIR = $(WORK_DIR)/utilities/DWCcode/error
|
|||
ARRAY_DIR = $(WORK_DIR)/utilities/DWCcode/arrays
|
||||
ASCII_DIR = $(WORK_DIR)/utilities/DWCcode/ascii
|
||||
STAT_DIR = $(WORK_DIR)/utilities/DWCcode/stat
|
||||
SPHERICAL_DIR = $(WORK_DIR)/utilities/DWCcode/spherical
|
||||
SORT_DIR = $(WORK_DIR)/utilities/DWCcode/sort
|
||||
SWITCH_DIR = $(WORK_DIR)/sbvar/switching
|
||||
VAR_DIR = $(WORK_DIR)/sbvar/var
|
||||
#################################################################################
|
||||
|
||||
# DW FILES
|
||||
INCLUDE_DIR := $(INCLUDE_DIR) -I$(MATRIX_DIR) -I$(ERROR_DIR) -I$(ARRAY_DIR) -I$(ASCII_DIR) -I$(STAT_DIR) -I$(SWITCH_DIR) -I$(VAR_DIR)
|
||||
VPATH := $(VPATH) $(MATRIX_DIR) $(ERROR_DIR) $(ARRAY_DIR) $(ASCII_DIR) $(STAT_DIR) $(SWITCH_DIR) $(VAR_DIR)
|
||||
OBJS := $(OBJS) bmatrix.o matrix.o dw_error.o dw_rand.o dw_matrix_rand.o dw_array.o dw_matrix_array.o dw_ascii.o dw_parse_cmd.o
|
||||
INCLUDE_DIR := -I$(MATRIX_DIR) -I$(ERROR_DIR) -I$(ARRAY_DIR) -I$(ASCII_DIR) -I$(STAT_DIR) -I$(SPHERICAL_DIR) -I$(SORT_DIR) -I$(SWITCH_DIR) -I$(VAR_DIR) $(INCLUDE_DIR)
|
||||
VPATH := $(VPATH) $(MATRIX_DIR) $(ERROR_DIR) $(ARRAY_DIR) $(ASCII_DIR) $(STAT_DIR) $(SPHERICAL_DIR) $(SORT_DIR) $(SWITCH_DIR) $(VAR_DIR)
|
||||
OBJS := $(OBJS) bmatrix.o matrix.o dw_error.o dw_rand.o dw_matrix_rand.o dw_array.o dw_matrix_array.o dw_matrix_sort.o dw_ascii.o dw_parse_cmd.o
|
||||
|
||||
# TAO FILES
|
||||
OBJS := $(OBJS)
|
||||
|
|
@ -189,29 +192,38 @@ OBJS := $(OBJS)
|
|||
# PROJECT FILE
|
||||
INCLUDE_DIR := $(INCLUDE_DIR)
|
||||
VPATH := $(VPATH)
|
||||
OBJS := $(OBJS) PrintDraws.o switch.o switchio.o VARbase.o VARio.o command_line_VAR.o
|
||||
OBJS_DRAWS := $(OBJS) PrintDraws.o switch.o switchio.o VARbase.o VARio.o command_line_VAR.o
|
||||
OBJS_TAO := tzmatlab.o mathlib.o cstz_dw.o
|
||||
OBJS_ESTIM := $(OBJS) $(OBJS_TAO) estimate.o VARbase.o VARio.o command_line_VAR.o switch.o switchio.o switch_opt.o csminwel.o
|
||||
OBJS_INIT := $(OBJS) create_init_file.o switch.o switchio.o VARbase.o VARio.o VARio_matlab.o
|
||||
OBJS_MHM_1 := $(OBJS) mhm_VAR_main_1.o mhm_VAR.o VARbase.o VARio.o command_line_VAR.o switch.o switchio.o
|
||||
OBJS_MHM_2 := $(OBJS) mhm_VAR_main_2.o spherical.o VARbase.o VARio.o switch.o switchio.o mhm_VAR.o
|
||||
OBJS_PROBA := $(OBJS) probabilities.o switch.o switchio.o VARbase.o VARio.o command_line_VAR.o
|
||||
|
||||
|
||||
|
||||
#################################################################################
|
||||
# OUTPUT
|
||||
|
||||
all: $(OUT_DIR)/sbvar_draws $(OUT_DIR)/sbvar_estimation $(OUT_DIR)/sbvar_init_file $(OUT_DIR)/sbvar_mhm_1 $(OUT_DIR)/sbvar_mhm_2 $(OUT_DIR)/sbvar_probabilities
|
||||
#all: $(OUT_DIR)/sbvar_draws $(OUT_DIR)/sbvar_estimation $(OUT_DIR)/sbvar_init_file $(OUT_DIR)/sbvar_mhm_1 $(OUT_DIR)/sbvar_mhm_2 $(OUT_DIR)/sbvar_probabilities
|
||||
all: $(OUT_DIR)/sbvar_draws $(OUT_DIR)/sbvar_init_file $(OUT_DIR)/sbvar_mhm_1 $(OUT_DIR)/sbvar_mhm_2 $(OUT_DIR)/sbvar_probabilities
|
||||
|
||||
$(OUT_DIR)/sbvar_draws: $(OBJS)
|
||||
$(OUT_DIR)/sbvar_draws: $(OBJS_DRAWS)
|
||||
$(CC) $(CFLAGS) $^ $(LIBS_DIR) $(LIBS) -o $(OUT_DIR)/sbvar_draws
|
||||
|
||||
$(OUT_DIR)/sbvar_estimation: $(OBJS)
|
||||
$(OUT_DIR)/sbvar_estimation: $(OBJS_ESTIM)
|
||||
$(CC) $(CFLAGS) $^ $(LIBS_DIR) $(LIBS) -o $(OUT_DIR)/sbvar_estimation
|
||||
|
||||
$(OUT_DIR)/sbvar_init_file: $(OBJS)
|
||||
$(OUT_DIR)/sbvar_init_file: $(OBJS_INIT)
|
||||
$(CC) $(CFLAGS) $^ $(LIBS_DIR) $(LIBS) -o $(OUT_DIR)/sbvar_init_file
|
||||
|
||||
$(OUT_DIR)/sbvar_mhm_1: $(OBJS)
|
||||
$(OUT_DIR)/sbvar_mhm_1: $(OBJS_MHM_1)
|
||||
$(CC) $(CFLAGS) $^ $(LIBS_DIR) $(LIBS) -o $(OUT_DIR)/sbvar_mhm_1
|
||||
|
||||
$(OUT_DIR)/sbvar_mhm_2: $(OBJS)
|
||||
$(OUT_DIR)/sbvar_mhm_2: $(OBJS_MHM_2)
|
||||
$(CC) $(CFLAGS) $^ $(LIBS_DIR) $(LIBS) -o $(OUT_DIR)/sbvar_mhm_2
|
||||
|
||||
$(OUT_DIR)/sbvar_probabilities: $(OBJS)
|
||||
$(OUT_DIR)/sbvar_probabilities: $(OBJS_PROBA)
|
||||
$(CC) $(CFLAGS) $^ $(LIBS_DIR) $(LIBS) -o $(OUT_DIR)/sbvar_probabilities
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,32 @@
|
|||
function [Ui,Vi,n0,np,ixmC0Pres] = upper_cholesky(nvar,nexo,options_ms)
|
||||
|
||||
lags = options_ms.nlags;
|
||||
indxC0Pres = options_ms.cross_restrictions;
|
||||
|
||||
Ui = cell(nvar,1);
|
||||
Vi = cell(nvar,1);
|
||||
n0 = zeros(nvar,1);
|
||||
np = zeros(nvar,1);
|
||||
if (nargin==2)
|
||||
nexo = 1;
|
||||
elseif (nargin==3)
|
||||
indxC0Pres = 0;
|
||||
end
|
||||
k = lags*nvar+nexo;
|
||||
Qi = zeros(nvar,nvar,nvar);
|
||||
Ri = zeros(k,k,nvar);
|
||||
|
||||
for ii=2:nvar
|
||||
Pi=diag(diag(ones(nvar-ii+1)),ii-1);
|
||||
Qi(:,:,ii-1)=Pi(1:nvar,1:nvar);
|
||||
Qi(:,:,nvar)=zeros(nvar,nvar);
|
||||
end
|
||||
|
||||
for n=1:nvar
|
||||
Ui{n} = null(Qi(:,:,n));
|
||||
Vi{n} = null(Ri(:,:,n));
|
||||
n0(n) = size(Ui{n},2);
|
||||
np(n) = size(Vi{n},2);
|
||||
end
|
||||
|
||||
ixmC0Pres = NaN;
|
||||
|
|
@ -1,22 +1,28 @@
|
|||
function [Ui,Vi,n0,np,ixmC0Pres] = upper_cholesky(lags,nvar,nexo,indxC0Pres)
|
||||
function [Ui,Vi,n0,np,ixmC0Pres] = lower_cholesky(nvar,nexo,options_ms)
|
||||
|
||||
lags = options_ms.nlags;
|
||||
indxC0Pres = options_ms.cross_restrictions;
|
||||
|
||||
Ui = cell(nvar,1);
|
||||
Vi = cell(nvar,1);
|
||||
n0 = zeros(nvar,1);
|
||||
np = zeros(nvar,1);
|
||||
|
||||
if (nargin==2)
|
||||
nexo = 1;
|
||||
nexo = 1;
|
||||
elseif (nargin==3)
|
||||
indxC0Pres = 0;
|
||||
end
|
||||
k = lags*nvar+nexo;
|
||||
|
||||
k = lags*nvar+nexo;
|
||||
Qi = zeros(nvar,nvar,nvar);
|
||||
Ri = zeros(k,k,nvar);
|
||||
|
||||
for ii=2:nvar
|
||||
Pi=diag(diag(ones(nvar-ii+1)),ii-1);
|
||||
Qi(:,:,ii-1)=Pi(1:nvar,1:nvar);
|
||||
Qi(:,:,nvar)=zeros(nvar,nvar);
|
||||
|
||||
Qi(ii-1,ii-1,ii)=1;
|
||||
Qi(:,:,ii)=Qi(:,:,ii)+Qi(:,:,ii-1);
|
||||
|
||||
end
|
||||
|
||||
for n=1:nvar
|
||||
|
|
|
|||
|
|
@ -71,7 +71,7 @@ IndxNmlr = [1 0 0 0 0 0]; % imported by nmlzvar.m
|
|||
%
|
||||
%options_.ms.indxgdls = 1; % 1: graph point forecast on variables; 0: disable
|
||||
nconstr1=nfqm; % number of the 1st set of constraints
|
||||
nconstr2=options_.ms.forecast ; % number of the 2nd set of constraints
|
||||
nconstr2=options_.forecast ; % number of the 2nd set of constraints
|
||||
nconstr=nconstr1+nconstr2; % q: 4 years -- 4*12 months.
|
||||
% When 0, no conditions directly on variables <<>>
|
||||
nconstr=0 %6*nconstr1;
|
||||
|
|
@ -274,7 +274,7 @@ if options_.ms.indxestima
|
|||
ye = [dateswd y];
|
||||
|
||||
%* Obtain linear restrictions
|
||||
[Uiconst,Viconst,n0,np,ixmC0Pres] = feval(options_.ms.restriction_fname,options_.ms.nlags ,nvar,nexo,options_.ms.cross_restrictions );
|
||||
[Uiconst,Viconst,n0,np,ixmC0Pres] = feval(options_.ms.restriction_fname,nvar,nexo,options_.ms );
|
||||
if min(n0)==0
|
||||
disp(' ')
|
||||
warning('A0: restrictions in dlrprior.m give no free parameter in one of equations')
|
||||
|
|
@ -574,7 +574,7 @@ if options_.ms.indxestima
|
|||
|
||||
%*** conditions in every period
|
||||
vpntM = fn_dataext(E1yrqm, E2yrqm,xdatae); % point value matrix with dates
|
||||
% vaveM = fn_dataext([yrEnd+1 0],[yrEnd+options_.ms.forecast 0],yact2yre); % average value matrix with dates
|
||||
% vaveM = fn_dataext([yrEnd+1 0],[yrEnd+options_.forecast 0],yact2yre); % average value matrix with dates
|
||||
for i=1:nconstr
|
||||
if i<=nconstr1
|
||||
valuecon(i) = vpntM(i,2+varcon(i)); % 2: first 2 elements are dates
|
||||
|
|
@ -593,10 +593,10 @@ if options_.ms.indxestima
|
|||
|
||||
% %*** average condtions over,say, options_.ms.freq periods.
|
||||
% if qmEnd==options_.ms.freq
|
||||
% vaveM = fn_dataext([yrEnd+1 0],[yrEnd+options_.ms.forecast 0],yact2yre); % average value matrix with dates
|
||||
% vaveM = fn_dataext([yrEnd+1 0],[yrEnd+options_.forecast 0],yact2yre); % average value matrix with dates
|
||||
% valuecon(1) = vaveM(1,2+varcon(1)); % 2: first 2 elements are dates
|
||||
% else
|
||||
% vaveM = fn_dataext([yrEnd 0],[yrEnd+options_.ms.forecast 0],yact2yre); % average value matrix with dates
|
||||
% vaveM = fn_dataext([yrEnd 0],[yrEnd+options_.forecast 0],yact2yre); % average value matrix with dates
|
||||
% yactrem = fn_dataext([yrEnd qmEnd+1],[yrEnd options_.ms.freq],xdatae);
|
||||
% valuecon(1) = sum(yactrem(:,2+varcon(1)),1)/length(stepcon{1});
|
||||
% % 2: first 2 elements are dates
|
||||
|
|
|
|||
|
|
@ -69,8 +69,8 @@ yrStart=options_.ms.initial_year;
|
|||
qmStart=options_.ms.initial_subperiod;
|
||||
yrEnd=options_.ms.final_year;
|
||||
qmEnd=options_.ms.final_subperiod;
|
||||
%options_.ms.forecast = 4; % number of years for forecasting
|
||||
if options_.ms.forecast<1
|
||||
%options_.forecast = 4; % number of years for forecasting
|
||||
if options_.forecast<1
|
||||
error('To be safe, the number of forecast years should be at least 1')
|
||||
end
|
||||
ystr=num2str(yrEnd);
|
||||
|
|
@ -82,7 +82,7 @@ if qmEnd==options_.ms.freq
|
|||
else
|
||||
E1yrqm = [yrEnd qmEnd+1]; % first year and quarter (month) after the sample
|
||||
end
|
||||
E2yrqm = [yrEnd+options_.ms.forecast qmEnd]; % end at the last month (quarter) of a calendar year after the sample
|
||||
E2yrqm = [yrEnd+options_.forecast qmEnd]; % end at the last month (quarter) of a calendar year after the sample
|
||||
[fdates,nfqm]=fn_calyrqm(options_.ms.freq,E1yrqm,E2yrqm); % forecast dates and number of forecast dates
|
||||
[sdates,nsqm] = fn_calyrqm(options_.ms.freq,[yrStart qmStart],[yrEnd qmEnd]);
|
||||
% sdates: dates for the whole sample (including options_.ms.nlags)
|
||||
|
|
@ -92,7 +92,7 @@ if nSample~=nsqm
|
|||
pause
|
||||
end
|
||||
imstp = 4*options_.ms.freq; % <<>> impulse responses (4 years)
|
||||
nayr = 4; %options_.ms.forecast; % number of years before forecasting for plotting.
|
||||
nayr = 4; %options_.forecast; % number of years before forecasting for plotting.
|
||||
|
||||
|
||||
%------- Prior, etc. -------
|
||||
|
|
@ -132,7 +132,7 @@ indxmsmdeqn = [0; 0; 0; 0]; %This option disenable using this in fn_rnrprior_co
|
|||
|
||||
|
||||
tdf = 3; % degrees of freedom for t-dist for initial draw of the MC loop
|
||||
nbuffer = 100; % a block or buffer of draws (buffer) that is saved to the disk (not memory)
|
||||
nbuffer = 1000; % a block or buffer of draws (buffer) that is saved to the disk (not memory)
|
||||
ndraws1=1*nbuffer; % 1st part of Monte Carlo draws
|
||||
ndraws2=10*ndraws1 % 2nd part of Monte Carlo draws
|
||||
seednumber = 0; %7910; %472534; % if 0, random state at each clock time
|
||||
|
|
|
|||
|
|
@ -10,14 +10,43 @@ addpath([swz_root '/mhm_specification']);
|
|||
|
||||
options.data = read_variables(options.datafile,options.varobs,[],options.xls_sheet,options.xls_range);
|
||||
|
||||
if options.forecast == 0
|
||||
options.forecast = 4;
|
||||
end
|
||||
|
||||
options.ms.output_file_tag = M.fname;
|
||||
%options.ms.markov_file = 'specification_2v2c.dat';
|
||||
%options.ms.mhm_file = 'MHM_input.dat';
|
||||
%options.ms.restriction_fname = 'ftd_upperchol3v';
|
||||
|
||||
if options.ms.upper_cholesky
|
||||
if options.ms.lower_cholesky
|
||||
error(['Upper Cholesky and lower Cholesky decomposition can''t be ' ...
|
||||
'requested at the same time!'])
|
||||
else
|
||||
options.ms.restriction_fname = 'upper_cholesky';
|
||||
end
|
||||
elseif options.ms.lower_cholesky
|
||||
options.ms.restriction_fname = 'lower_cholesky';
|
||||
elseif ~isempty(options.ms.Qi) && ~isempty(options.ms.Ri)
|
||||
options.ms.restriction_fname = 'exclusions';
|
||||
end
|
||||
|
||||
if ms_flag == 1
|
||||
sz_prd(options)
|
||||
% changing some option names to match SWZ code
|
||||
options.ms.firstMetrop = options.ms.draws_nbr_burn_in_1;
|
||||
options.ms.secondMetrop = options.ms.draws_nbr_burn_in_2;
|
||||
options.ms.ndrawsmv = options.ms.draws_nbr_mean_var_estimate;
|
||||
options.ms.ndrawsmhm = options.ms.draws_nbr_modified_harmonic_mean;
|
||||
options.ms.tfmhm = options.ms.thinning_factor;
|
||||
options.ms.svd = options.ms.dirichlet_scale;
|
||||
% are these options necessary ?
|
||||
options.ms.opt1 = 1;
|
||||
options.ms.opt2 = 1;
|
||||
options.ms.opt3 = 1;
|
||||
% temporary fix
|
||||
options.ms.markov_file = 'markov_file';
|
||||
sz_prd(M,options);
|
||||
else
|
||||
swz_mardd(options);
|
||||
end
|
||||
|
|
|
|||
|
|
@ -0,0 +1,126 @@
|
|||
function swz_write_markov_file(fname,M,options)
|
||||
|
||||
n_chains = length(options.ms.ms_chain);
|
||||
nvars = size(options.varobs,1);
|
||||
|
||||
fh = fopen(fname,'w');
|
||||
%/******************************************************************************/
|
||||
%/********************* Markov State Variable Information **********************/
|
||||
%/******************************************************************************/
|
||||
|
||||
fprintf(fh,'//== Flat Independent Markov States and Simple Restrictions ==//\n\n');
|
||||
|
||||
|
||||
%//This number is NOT used but read in.
|
||||
fprintf(fh,'//== Number Observations ==//\n');
|
||||
fprintf(fh,'0\n\n');
|
||||
|
||||
fprintf(fh,'//== Number Independent State Variables ==//\n');
|
||||
fprintf(fh,'%d\n\n',n_chains);
|
||||
|
||||
for i_chain = 1:n_chains
|
||||
|
||||
%//=====================================================//
|
||||
%//== state_variable[i] (1 <= i <= n_state_variables) ==//
|
||||
%//=====================================================//
|
||||
fprintf(fh,'//== Number of states for state_variable[%d] ==//\n', ...
|
||||
i_chain);
|
||||
n_states = length(options.ms.ms_chain(i_chain).state);
|
||||
fprintf(fh,'%d\n\n',n_states);
|
||||
|
||||
%//== 03/15/06: DW TVBVAR code reads the data below and overwrite the prior data read somewhere else if any.
|
||||
%//== Each column contains the parameters for a Dirichlet prior on the corresponding
|
||||
%//== column of the transition matrix. Each element must be positive. For each column,
|
||||
%//== the relative size of the prior elements determine the relative size of the elements
|
||||
%//== of the transition matrix and overall larger sizes implies a tighter prior.
|
||||
fprintf(fh,['//== Transition matrix prior for state_variable[%d]. ' ...
|
||||
'(n_states x n_states) ==//\n'],i_chain);
|
||||
Alpha = ones(n_states,n_states);
|
||||
for i_state = 1:n_states
|
||||
p = 1-1/options.ms.ms_chain(i_chain).state(i_state).duration;
|
||||
Alpha(i_state,i_state) = p*(n_states-1)/(1-p);
|
||||
fprintf(fh,'%22.16f',Alpha(i_state,:));
|
||||
fprintf(fh,'\n');
|
||||
end
|
||||
|
||||
fprintf(fh,['\n//== Free Dirichet dimensions for state_variable[%d] ' ...
|
||||
'==//\n'],i_chain);
|
||||
fprintf(fh,'%d ',repmat(n_states,1,n_states));
|
||||
fprintf(fh,'\n\n');
|
||||
|
||||
%//== The jth restriction matrix is n_states-by-free[j]. Each row of the restriction
|
||||
%//== matrix has exactly one non-zero entry and the sum of each column must be one.
|
||||
fprintf(fh,['//== Column restrictions for state_variable[%d] ' ...
|
||||
'==//\n'],i_chain);
|
||||
if n_states == 2
|
||||
M = eye(2);
|
||||
for i_state = 1:2
|
||||
for j_state = 1:2
|
||||
fprintf(fh,'%d ',M(j_state,:));
|
||||
fprintf(fh,'\n');
|
||||
end
|
||||
fprintf(fh,'\n');
|
||||
end
|
||||
else
|
||||
end
|
||||
end
|
||||
|
||||
%/******************************************************************************/
|
||||
%/******************************* VAR Parameters *******************************/
|
||||
%/******************************************************************************/
|
||||
%//NOT read
|
||||
fprintf(fh,'//== Number Variables ==//\n');
|
||||
fprintf(fh,'%d\n\n',nvars);
|
||||
|
||||
%//NOT read
|
||||
fprintf(fh,'//== Number Lags ==//\n');
|
||||
fprintf(fh,'%d\n\n',options.ms.nlags);
|
||||
|
||||
%//NOT read
|
||||
fprintf(fh,'//== Exogenous Variables ==//\n');
|
||||
fprintf(fh,'1\n\n');
|
||||
|
||||
|
||||
%//== nvar x n_state_variables matrix. In the jth row, a non-zero value implies that
|
||||
%this state variable controls the jth column of A0 and Aplus
|
||||
fprintf(fh,['//== Controlling states variables for coefficients ==//\' ...
|
||||
'n']);
|
||||
|
||||
for i_var = 1:nvars
|
||||
for i_chain = 1:n_chains
|
||||
if ~isfield(options.ms.ms_chain(i_chain),'svar_coefficients')
|
||||
i_equations = 0;
|
||||
else
|
||||
i_equations = ...
|
||||
options.ms.ms_chain(i_chain).svar_coefficients.equations;
|
||||
end
|
||||
if strcmp(i_equations,'ALL') || any(i_equations == i_var)
|
||||
fprintf(fh,'%d ',1);
|
||||
else
|
||||
fprintf(fh,'%d ',0);
|
||||
end
|
||||
end
|
||||
fprintf(fh,'\n');
|
||||
end
|
||||
|
||||
%//== nvar x n_state_variables matrix. In the jth row, a non-zero value implies that
|
||||
%this state variable controls the jth diagonal element of Xi
|
||||
fprintf(fh,'\n//== Controlling states variables for variance ==//\n');
|
||||
for i_var = 1:nvars
|
||||
for i_chain = 1:n_chains
|
||||
if ~isfield(options.ms.ms_chain(i_chain),'svar_variances')
|
||||
i_equations = 0;
|
||||
else
|
||||
i_equations = ...
|
||||
options.ms.ms_chain(i_chain).svar_variances.equations;
|
||||
end
|
||||
if strcmp(i_equations,'ALL') || any(i_equations == i_var)
|
||||
fprintf(fh,'%d ',1);
|
||||
else
|
||||
fprintf(fh,'%d ',0);
|
||||
end
|
||||
end
|
||||
fprintf(fh,'\n');
|
||||
end
|
||||
|
||||
fclose(fh);
|
||||
|
|
@ -0,0 +1,38 @@
|
|||
function swz_write_mhm_input(fname,options_ms)
|
||||
|
||||
fh = fopen(fname,'w');
|
||||
|
||||
|
||||
fprintf(fh,'/**********************************************************\n');
|
||||
fprintf(fh,' *** This input file is read by swzmsbvar_mhm_1 and swzmsbvar_mhm_1.exe only, NOT by swzmsbvar_printdraws.exe.\n');
|
||||
fprintf(fh,' ***\n');
|
||||
fprintf(fh,' **********************************************************/\n');
|
||||
|
||||
fprintf(fh,'\n\n//------------- 1st set of posterior draws to find optimal scales for Metropolis (30000). ---------------\n');
|
||||
fprintf(fh,'//== number draws for first burn-in ==// //For determining the Metropolis scales only.\n');
|
||||
fprintf(fh,'%d\n\n',options_ms.draws_nbr_burn_in_1);
|
||||
|
||||
fprintf(fh,'//------------- MCMC burn-in draws once the Metropolis scales (previous stage) are fixed. --------------\n');
|
||||
fprintf(fh,'//------------- 2nd set of standard burn-in posterior draws to throw away the initial draws (10000). ---------------\n');
|
||||
fprintf(fh,'//== number draws for second burn-in ==//\n');
|
||||
fprintf(fh,'%d\n\n',options_ms.draws_nbr_burn_in_2);
|
||||
|
||||
fprintf(fh,'//--------------- 1st set of posterior draws to compute the mean and variance for the weighting function in the MHM (200000) ----------------\n');
|
||||
fprintf(fh,'//== number draws to estimate mean and variance ==//\n');
|
||||
fprintf(fh,'%d\n\n',options_ms.draws_nbr_mean_var_estimate);
|
||||
|
||||
fprintf(fh,'//--------------- Only applied to mhm_2 process: total number of MCMC draws = thinning factor * 2nd set of saved posterior draws ----------------\n');
|
||||
fprintf(fh,'//== thinning factor for modified harmonic mean process ==//\n');
|
||||
fprintf(fh,'%d\n\n',options_ms.thinning_factor);
|
||||
|
||||
fprintf(fh,'//--------------- 2nd set of saved posterior draws from MHM_2 (second stage): saved draws AFTER thinning (1000000) ----------------\n');
|
||||
fprintf(fh,'//== number draws for modified harmonic mean process ==//\n');
|
||||
fprintf(fh,'%d\n\n',options_ms.draws_nbr_modified_harmonic_mean);
|
||||
|
||||
fprintf(fh,'//------- 1st stage: computing all three tightness factors for Dirichlet. ---------\n');
|
||||
fprintf(fh,'//------- 2nd stage: hard-code the second scale factor (in principle, we can do all three). ---------\n');
|
||||
fprintf(fh,'//------- It seems that Dan''s code only use the first element of the following scales. The scale applies to the Dirichlet''s hyperparameter alpha for the diagonal of the transition matrix in the weighting function. Note that the weighting function for the transition matrix parameters is Dirichlet. ---------\n');
|
||||
|
||||
fprintf(fh,'//== scale values for Dirichlet distribution ==//\n');
|
||||
fprintf(fh,'1\n\n');
|
||||
fprintf(fh,'%d\n',options_ms.dirichlet_scale);
|
||||
|
|
@ -1,4 +1,4 @@
|
|||
function sz_prd(options_)
|
||||
function sz_prd(M,options_)
|
||||
%==========================================================================
|
||||
%== Directory structure
|
||||
%==========================================================================
|
||||
|
|
@ -53,7 +53,7 @@ markov_file = [options_.ms.markov_file '.dat'];
|
|||
%options_.ms.mhm_file = 'MHM_input.dat';
|
||||
|
||||
|
||||
mhm_file = [mhm_spec_path '/MHM_input.dat'];
|
||||
mhm_file = '/MHM_input.dat';
|
||||
%options_.ms.proposal_draws = 100000;
|
||||
|
||||
%==========================================================================
|
||||
|
|
@ -277,7 +277,7 @@ if options_.ms.create_initialization_file == 1
|
|||
Ui{kj} = eye(nvar); Vi{kj} = eye(ncoef);
|
||||
end
|
||||
else
|
||||
eval(['[Ui,Vi,n0,np,ixmC0Pres] = ' options_.ms.restriction_fname '(options_.ms.nlags,nvar,nexo,indxC0Pres);'])
|
||||
[Ui,Vi,n0,np,ixmC0Pres] = feval(options_.ms.restriction_fname,nvar,nexo,options_.ms);
|
||||
if min(n0)==0
|
||||
disp('A0: restrictions give no free parameters in one of equations')
|
||||
return
|
||||
|
|
@ -326,7 +326,7 @@ if options_.ms.create_initialization_file == 1
|
|||
%
|
||||
tic
|
||||
[fhat,xhat,grad,Hhat,itct,fcount,retcodehat] = csminwel('fn_a0freefun',x,H0,'fn_a0freegrad',crit,nit,Ui,nvar,n0,fss,H0inv);
|
||||
endtime = toc
|
||||
endtime = toc;
|
||||
|
||||
A0hat = fn_tran_b2a(xhat,Ui,nvar,n0);
|
||||
|
||||
|
|
@ -480,9 +480,11 @@ if options_.ms.create_initialization_file == 1
|
|||
%======================================================================
|
||||
%== Create C initialization filename
|
||||
%======================================================================
|
||||
swz_write_markov_file(markov_file,M,options_)
|
||||
if use_linux == 1
|
||||
create_init_file=[c_path,'/sbvar_init_file ',matlab_filename,' ',m_spec_path,'/',markov_file,' ',options_.ms.output_file_tag];
|
||||
system(create_init_file) %Run operating system command and return result
|
||||
create_init_file=[c_path,'/sbvar_init_file ',matlab_filename,' ',markov_file,' ',options_.ms.output_file_tag];
|
||||
system(create_init_file); %Run operating system command and return
|
||||
%result
|
||||
else
|
||||
create_init_file=[c_path,'\sbvar_init.exe ',matlab_filename,' ',m_spec_path,'\',markov_file,' ',options_.ms.output_file_tag];
|
||||
dos(create_init_file)
|
||||
|
|
@ -495,10 +497,10 @@ end
|
|||
%==========================================================================
|
||||
if options_.ms.estimate_msmodel == 1
|
||||
if use_linux == 1
|
||||
perform_estimation=[c_path,'/sbvar_estimate -ft ',options_.ms.output_file_tag];
|
||||
system(perform_estimation)
|
||||
perform_estimation=[c_path,'/sbvar_estimation -ft ',options_.ms.output_file_tag];
|
||||
% system(perform_estimation);
|
||||
else
|
||||
perform_estimation=[c_path,'\sbvar_estimate.exe -ft ',options_.ms.output_file_tag];
|
||||
perform_estimation=[c_path,'\sbvar_estimation.exe -ft ',options_.ms.output_file_tag];
|
||||
dos(perform_estimation)
|
||||
end
|
||||
end
|
||||
|
|
@ -508,13 +510,15 @@ end
|
|||
%== Compute marginal data density
|
||||
%==========================================================================
|
||||
if options_.ms.compute_mdd == 1
|
||||
mhm_file = ['mhm_input_' M.fname '.dat'];
|
||||
swz_write_mhm_input(mhm_file,options_.ms);
|
||||
if use_linux == 1
|
||||
compute_mdd1=[c_path,'/sbvar_mhm_1 -ft ',options_.ms.output_file_tag,' -fi ',mhm_spec_path,'/',mhm_file];
|
||||
compute_mdd1=[c_path,'/sbvar_mhm_1 -ft ',options_.ms.output_file_tag,' -fi ',mhm_file];
|
||||
system(compute_mdd1);
|
||||
compute_mdd2=[c_path,'/sbvar_mhm_2 -ft ',options_.ms.output_file_tag,' -d ',int2str(options_.ms.proposal_draws),' -t 3'];
|
||||
system(compute_mdd2);
|
||||
else
|
||||
compute_mdd1=[c_path,'\sbvar_mhm_1.exe -ft ',options_.ms.output_file_tag,' -fi ',mhm_spec_path,'/',mhm_file];
|
||||
compute_mdd1=[c_path,'\sbvar_mhm_1.exe -ft ',options_.ms.output_file_tag,' -fi ',mhm_file];
|
||||
system(compute_mdd1);
|
||||
compute_mdd2=[c_path,'\sbvar_mhm_2.exe -ft ',options_.ms.output_file_tag,' -d ',int2str(options_.ms.proposal_draws),' -t 3'];
|
||||
system(compute_mdd2);
|
||||
|
|
|
|||
|
|
@ -1100,19 +1100,20 @@ SvarIdentificationStatement::writeOutput(ostream &output, const string &basename
|
|||
if (!upper_cholesky_present && !lower_cholesky_present)
|
||||
{
|
||||
int n = symbol_table.endo_nbr();
|
||||
int m = symbol_table.exo_nbr();
|
||||
int r = getMaxLag();
|
||||
// int m = symbol_table.exo_nbr();
|
||||
int m = 1; // this is the constant, not the shocks
|
||||
int r = getMaxLag()+1;
|
||||
int k = r*n+m;
|
||||
|
||||
if (k < 1)
|
||||
{
|
||||
cerr << "ERROR: lag = " << r
|
||||
<< ", number of endogenous variables = " << n
|
||||
<< ", number of exogenous variables = " << m
|
||||
<< ". If this is not a logical error in the specification"
|
||||
<< " of the .mod file, please report it to the Dynare Team." << endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
{
|
||||
cerr << "ERROR: lag = " << r
|
||||
<< ", number of endogenous variables = " << n
|
||||
<< ", number of exogenous variables = " << m
|
||||
<< ". If this is not a logical error in the specification"
|
||||
<< " of the .mod file, please report it to the Dynare Team." << endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
if (n < 1)
|
||||
{
|
||||
cerr << "ERROR: Number of endogenous variables = " << n << "< 1. If this is not a logical "
|
||||
|
|
@ -1140,7 +1141,7 @@ SvarIdentificationStatement::writeOutput(ostream &output, const string &basename
|
|||
}
|
||||
|
||||
if (it->first.first == 0)
|
||||
output << "options_.ms.Qi(" << h+1 << ", " << j << ", "<< i << ");" << endl;
|
||||
output << "options_.ms.Qi(" << h+1 << ", " << j << ", "<< i << ") = 1;" << endl;
|
||||
else if (it->first.first > 0)
|
||||
{
|
||||
if ((it->first.first-1)*n+j > k)
|
||||
|
|
@ -1148,7 +1149,7 @@ SvarIdentificationStatement::writeOutput(ostream &output, const string &basename
|
|||
cerr << "ERROR: lag =" << it->first.first << ", num endog vars = " << n << "current endog var index = " << j << ". Index "
|
||||
<< "out of bounds. If the above does not represent a logical error, please report this to the Dyanre Team." << endl;
|
||||
}
|
||||
output << "options_.ms.Ri(" << h+1 << ", " << (it->first.first-1)*n+j << ", "<< i << ");" << endl;
|
||||
output << "options_.ms.Ri(" << h+1 << ", " << (it->first.first-1)*n+j << ", "<< i << ") = 1;" << endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
|
|
|
|||
|
|
@ -0,0 +1,16 @@
|
|||
addpath '/home/michel/dynare/svn/dynare/trunk/matlab/swz';
|
||||
var Y Pie R;
|
||||
|
||||
model;
|
||||
Y = 0;
|
||||
Pie = 0;
|
||||
R = 0;
|
||||
end;
|
||||
|
||||
varobs Y Pie R;
|
||||
|
||||
|
||||
sbvar(datafile = data,freq=4,initial_year=1959,final_year=2005,nlags=4,restriction_fname=upper_cholesky);//for SBVAR
|
||||
|
||||
//ms_sbvar(datafile = data,freq=4,initial_year=1959,final_year=2005,nlags=4,restriction_fname=ftd_upperchol3v,
|
||||
// markov_file=specification_2v2c,mhm_file=MHM_input); //for markov switching
|
||||
Loading…
Reference in New Issue