Adding Mex Functionality to work with new MS-SBVAR Code, includes...

Adding Mex Functionality to work with new MS-SBVAR Code, includes Forecasting,IRF,Variance Decomposition, Plotting and new global settings

matlab/global_initialization.m

@@ -255,66 +255,70 @@ options_.initial_date.subperiod = 0;
 % discretionary policy
 options_.discretionary_policy = 0;
 
-% SWZ SBVAR
+%% SWZ SBVAR
+options_.ms.filtered_probabilities = 1;
+options_.ms.real_time_smoothed_probabilities = 0;
 options_.ms.freq = 1;
 options_.ms.initial_subperiod = 1;
 options_.ms.final_subperiod=4;
-options_.ms.log_var = [ ];
-options_.ms.forecast = 1;
 options_.ms.nlags = 1;
+options_.ms.restriction_fname = 0;
 options_.ms.cross_restrictions = 0;
 options_.ms.contemp_reduced_form = 0;
-options_.ms.real_pseudo_forecast = 0;
 options_.ms.bayesian_prior = 1;
-options_.ms.dummy_obs = 0;
-options_.ms.ncsk = 0;
-options_.ms.nstd = 6;
-options_.ms.ninv = 1000;
-options_.ms.indxparr = 1;
-options_.ms.indxovr = 0;
-options_.ms.aband = 1;
-options_.ms.indxap = 1;
-options_.ms.apband = 1;
-options_.ms.indximf = 1;
-options_.ms.imfband = 1;
-options_.ms.indxfore = 0;
-options_.ms.foreband = 0;
-options_.ms.indxgforhat = 1;
-options_.ms.indxgimfhat = 1;
-options_.ms.indxestima = 1;
-options_.ms.indxgdls = 1;
-options_.ms.cms =0;
-options_.ms.ncms = 0;
-options_.ms.eq_cms = 1;
-options_.ms.cnum = 0;
-options_.ms.banact = 1;
-options_.ms.nstates = 2;
-options_.ms.indxscalesstates = 0;
 options_.ms.alpha = 1.0;
 options_.ms.beta = 1.0;
 options_.ms.gsig2_lmd = 50^2;
 options_.ms.gsig2_lmdm = 50^2;
-options_.ms.q_diag = 0.85;
-options_.ms.flat_prior = 0;
+options_.ms.lower_cholesky = 0;
+options_.ms.upper_cholesky = 0;
 options_.ms.create_initialization_file = 1;
-options_.ms.estimate_msmodel = 1;
 options_.ms.compute_mdd = 1;
 options_.ms.compute_probabilities = 1;
-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.coefficients_prior_hyperparameters = [1.0 1.0 0.1 1.2 1.0 1.0];
+options_.ms.irf = 0;
+options_.ms.bayesian_irf = 0;
+options_.ms.forecast = 0;
+options_.ms.variance_decomposition = 0; 
 options_.ms.thinning_factor = 1;
-options_.ms.dirichlet_scale = [1.0 1.5 2.0];
-options_.ms.vlistlog = [];
+options_.ms.shock_draws = 10000;
+options_.ms.shocks_per_parameter = 10;
+options_.ms.percentiles = [.16 .5 .84];
+options_.ms.mode_compute = 1;
+options_.ms.mode_file = 0;
+options_.ms.estimate.random = 5;
+options_.ms.estimate.random_max = 20;
+options_.ms.estimate.random_tol_obj = 0.1;
+options_.ms.estimate.random_tol_parms = 0.1;
+options_.ms.estimate.cb = 1e-3;
+options_.ms.estimate.ce = 1e-6;
+options_.ms.estimate.ci = 0.1;
+options_.ms.estimate.ib = 50;
+options_.ms.estimate.ii = 2.0;
+options_.ms.estimate.mb = 100;
+options_.ms.load_mh_file = 0;
+options_.ms.mh_replic = 20000;
+options_.ms.thinning_factor = 10;
+options_.ms.drop = 2000;
+options_.ms.adapative_mh_draws = 30000;
+options_.ms.mdd_proposal_type = [3 0.1 0.9];
+options_.ms.mdd_use_mean_center = 0;
+options_.ms.mdd_proposal_draws = 100000;
+options_.ms.filtered_probabilities = 0;
+options_.ms.real_time_smoothed_probabilities = 0;
+options_.ms.irf_shocks_per_parameter = options_.ms.shocks_per_parameter;
+options_.ms.irf_shock_draws = options_.ms.shock_draws;
+options_.ms.irf_thinning_factor = 1;
+options_.ms.irf_filtered = 0;
+options_.ms.forecast_shocks_per_parameter = options_.ms.shocks_per_parameter;
+options_.ms.forecast_shock_draws = options_.ms.shock_draws;
+options_.ms.forecast_thinning_factor = 1;
+options_.ms.forecast_data_obs = 0;
+options_.ms.vd_shocks_per_parameter = options_.ms.shocks_per_parameter;
+options_.ms.vd_shock_draws = options_.ms.shock_draws;
+options_.ms.vd_thinning_factor = 1;
+options_.ms.vd_filtered = 0;
+
 
 % Shock decomposition
 options_.parameter_set = [];

matlab/ms-sbvar/ms_sbvar.m

@@ -35,17 +35,6 @@ elseif ~isempty(options.ms.Qi) && ~isempty(options.ms.Ri)
 end
 
 if ms_flag == 1
-    % changing some option names to match MS 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);

matlab/ms-sbvar/ms_write_markov_file.m

@@ -15,7 +15,7 @@ function ms_write_markov_file(fname,M,options)
     fprintf(fh,'//== Number Observations ==//\n');
     fprintf(fh,'0\n\n');
 
-    fprintf(fh,'//== Number Independent State Variables ==//\n');
+    fprintf(fh,'//== Number independent state variables ==//\n');
     fprintf(fh,'%d\n\n',n_chains);
 
     for i_chain = 1:n_chains

matlab/ms-sbvar/plot_ms_forecast.m

+function plot_ms_forecast(forecast,title_)
+% function [] = plot_ms_forecast(forecast,names)
+% plots the forecast from the output from a ms-sbvar
+%
+% INPUTS
+%   forecast should be in the form (percentile x horizon x nvar ), if banded otherwise
+%     ( horizon x nvar )
+%
+%   title: optional super title
+%
+
+% Copyright (C) 2011 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 M_ oo_ options_
+
+  nc = 2;
+  nr = 2;
+  nvars = M_.endo_nbr;
+  endo_names = M_.endo_names;
+  fname = M_.fname;
+  
+  var_list = endo_names(1:M_.orig_endo_nbr,:);
+
+  i_var = [];
+  for i = 1:size(var_list)
+      tmp = strmatch(var_list(i,:),endo_names,'exact');
+      if isempty(tmp)
+          error([var_list(i,:) ' isn''t and endogenous variable'])
+      end
+      i_var = [i_var; tmp];
+  end
+  nvar = length(i_var);
+  
+  dims = size(forecast);
+  
+  if nargin < 2
+    title_ = '';
+  end
+  
+  if (length(dims) == 2)
+    % Point Forecast (horizon x nvars )
+    horizon = dims(1);
+    num_percentiles = 1;
+  elseif (length(dims) == 3)
+    % Banded Forecast
+    horizon = dims(2);
+    num_percentiles = dims(1);  
+  else
+    error('The impulse response matrix passed to be plotted does not appear to be the correct size');
+  end
+
+  if ~exist(fname, 'dir')
+      mkdir('.',fname);
+  end
+  if ~exist([fname '/graphs'])
+      mkdir(fname,'graphs');
+  end
+
+  n_fig = 1;
+  if num_percentiles == 1
+    plot_point_forecast(forecast,nvars,endo_names);
+    
+  else
+    plot_banded_forecast(forecast,nvars,endo_names,num_percentiles);
+  end
+
+  function plot_point_forecast(forecast,nvars,names)
+    fig = figure('Name','Forecast (I)'); 
+    m = 1;
+    for j=1:nvars
+      if m > nr*nc
+        %eval(['print -deps ' fname '/graphs/forcst' int2str(n_fig) '.eps;'])
+        n_fig =n_fig+1;
+        eval(['figure(''Name'',''Forecast (' int2str(n_fig) ')'');']);
+        m = 1;
+      end
+      subplot(nr,nc,m);
+      vn = deblank(names(i_var(j),:));
+      plot(forecast(:,j))
+      title(vn,'Interpreter','none');
+      %suptitle(title_);
+      grid on;
+      m = m+1;
+    end
+    if m > 1
+        %eval(['print -deps ' fname '/graphs/forcst' int2str(n_fig) '.eps;'])
+    end
+  end
+  
+  function plot_banded_forecast(forecast,nvars,names,num_percentiles)
+    fig = figure('Name',[title_ ' (1)']); 
+    %suptitle(title_);
+    m = 1;
+    
+    for j=1:nvars
+      if m > nr*nc
+        suptitle(title_);
+        %eval(['print -deps ' fname '/graphs/forcst' int2str(n_fig) '.eps;'])
+        n_fig =n_fig+1;
+        eval(['figure(''Name'',''' title_ ' (' int2str(n_fig) ')'');']);
+        m = 1;
+      end
+      subplot(nr,nc,m);
+      vn = deblank(names(i_var(j),:));
+      for k=1:num_percentiles
+        if ceil(num_percentiles/2) == k
+          plot(forecast(k,:,j),'LineWidth',1.5)
+        else
+          plot(forecast(k,:,j),'LineWidth',1.1)
+        end
+        if k==1
+          hold on;
+        end
+      end
+      title(vn,'Interpreter','none');
+      hold off
+      grid on;
+      m = m+1;
+    end
+    if m > 1
+        %eval(['print -deps ' fname '/graphs/forcst' int2str(n_fig) '.eps;'])
+    end
+  end
+  
+
+end
\ No newline at end of file

matlab/ms-sbvar/plot_ms_irf.m

+function plot_ms_irf(irf,names,title_)
+% function [] = plot_ms_irf(irf,names)
+% plots the impulse responses from the output from a ms-sbvar
+%
+% INPUTS
+%   irf should be in the form (percentile x horizon x (nvar x nvar)), if banded otherwise
+%     ( horizon x (nvar x nvar) )
+%
+%   names: character list of the names of the variables
+%
+%   title: optional super title
+%
+% The element in position (k,i+j*nvars) of ir is the response of the ith 
+% variable to the jth shock at horizon k.  Horizon 0 is the contemporaneous 
+% response.
+
+% Copyright (C) 2011 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 M_ oo_ options_
+
+  nc = 2;
+  nr = 2;
+  nvars = M_.endo_nbr;
+  endo_names = M_.endo_names;
+  fname = M_.fname;
+  
+  var_list = endo_names(1:M_.orig_endo_nbr,:);
+
+  i_var = [];
+  for i = 1:size(var_list)
+      tmp = strmatch(var_list(i,:),endo_names,'exact');
+      if isempty(tmp)
+          error([var_list(i,:) ' isn''t and endogenous variable'])
+      end
+      i_var = [i_var; tmp];
+  end
+
+  nvar = length(i_var);
+
+  dims = size(irf);
+  
+  if nargin < 3
+    title_ = '';
+  end
+  
+  if (length(dims) == 2)
+    % Point IRF (horizon x (nvarsxnvars) )
+    horizon = dims(1);
+    num_percentiles = 1;
+  elseif (length(dims) == 3)
+    % Banded IRF
+    horizon = dims(2);
+    num_percentiles = dims(1);
+  else
+    error('The impulse response matrix passed to be plotted does not appear to be the correct size');
+  end
+  
+  if size(names,1) ~= nvars
+    error('The names passed are not the same length as the number of variables')
+  end
+  
+  if ~exist(fname, 'dir')
+      mkdir('.',fname);
+  end
+  if ~exist([fname '/graphs'])
+      mkdir(fname,'graphs');
+  end
+  
+  
+  if num_percentiles == 1
+    % loop through the shocks
+    for s=1:nvars
+      shock = zeros(horizon,nvars);
+      for i=1:nvars
+        shock(:,i) = irf(:,((i-1) + ((s-1)*nvars)+1));
+      end
+      plot_point_irf_for_shock(shock,nvars,names,names(s,:));
+    end
+  else
+    for s=1:nvars
+      shock = zeros(horizon,nvars,num_percentiles);
+      for n=1:num_percentiles
+        for i=1:nvars
+          shock(:,i,n) = irf(n,:,((i-1) + ((s-1)*nvars)+1));
+        end
+      end
+      plot_banded_irf_for_shock(shock,nvars,names,names(s,:));
+    end
+  end
+  
+  function [fig] = plot_point_irf_for_shock(irf,nvars,names,shock_name)
+    fig = figure('Name',title_);
+    for k=1:nvars
+      subplot(ceil(sqrt(nvars)), ceil(sqrt(nvars)),k);
+      plot(irf(:,k))
+      title([names(k,:) ' shock from ' shock_name]);
+    end
+    %suptitle(title_);
+  end
+  
+  function [fig] = plot_banded_irf_for_shock(irf,nvars, names, shock_name)
+    fig = figure('Name',title_);
+    npercentiles = size(irf,3);
+    for ii=1:nvars
+      subplot(ceil(sqrt(nvars)), ceil(sqrt(nvars)),ii);
+      for nn=1:npercentiles
+        plot(irf(:,ii,nn))
+        hold on
+      end
+      hold off
+      title([names(ii,:) ' shock from ' shock_name]);
+    end
+    %suptitle(title_);
+  end
+  
+end
\ No newline at end of file

matlab/ms-sbvar/plot_ms_probabilities.m

+function plot_ms_probabilities(computed_probabilities)
+% function plot_ms_probabilities(computed_probabilities)
+% Plots the regime probablities for each graph
+%
+% INPUTS
+%    computed_probabilities:      Txnstates 
+%
+% OUTPUTS
+%    none
+%
+% SPECIAL REQUIREMENTS
+%    none
+
+% Copyright (C) 2011 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 M_ oo_ options_
+  
+  [T,num_grand_regimes] = size(computed_probabilities);
+  num_chains = length(options_.ms.ms_chain);
+  for i=1:num_chains
+    chains(i).num_regimes = length(options_.ms.ms_chain(i).state);
+    chains(i).probabilities = zeros([T,chains(i).num_regimes]);
+  end
+  
+  for t=1:T
+    chains = iterate_chain(computed_probabilities(t,:), t, chains, 1, num_chains);
+  end
+  
+  for i=1:num_chains
+    figure
+    plot(chains(i).probabilities,'LineWidth', 1.2);
+    ltxt = {};
+    for j=1:chains(i).num_regimes
+      ltxt{j} = ['Regime ' num2str(j)];
+    end
+    legend(ltxt{:});
+    title(['Chain ' num2str(i)]);
+    ylim([0 1.0]);
+  end
+  
+  
+end
+
+function [chains] = iterate_chain(probs, t, chains, chain, num_chains)
+  offset_length = length(probs)/chains(chain).num_regimes;
+  for i=1:chains(chain).num_regimes
+    p = probs( (i-1)*offset_length+1 : i*offset_length );
+    chains(chain).probabilities(t, i) = chains(chain).probabilities(t, i) + sum( p );
+    if chain < num_chains
+      chains = iterate_chain(p, t, chains, chain+1, num_chains);
+    end
+  end
+end
\ No newline at end of file

matlab/ms-sbvar/plot_ms_variance_decomposition.m

+function plot_ms_variance_decomposition(vd, title_, varargin)
+% plot the variance decomposition of shocks
+%
+% Inputs
+%		shocks: matrix of the individual shocks Tx(KxK)with J=number of shocks
+%
+% Optional Inputs
+%		'data': the actual data, TxK with K=number of data series
+%		'steady': the steady state value, TxK
+%		'shock_names': to specify the names of the shocks
+%		'series_names': to specify the names of the different series
+%		'dates': pass a date vector to use, otherwise will just index on 1:T
+%		'colors': Jx3 list of the rgb colors to use for each shock
+%
+% Example:
+% plot_historic_decomposition(shocks,'VD','shock_names',shock_names,'series_names',series_names)
+
+% Copyright (C) 2011 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 M_ oo_ options_
+  nvars = M_.endo_nbr;
+  endo_names = M_.endo_names;
+
+	dims = size(vd);
+	if length(dims) == 3
+  	[T,K,J] = dims;
+  	shocks = vd;
+  else
+    T = dims(1);
+    K = nvars;
+    J = nvars;
+    temp_vd = zeros(T,K,J);
+    for i=1:nvars
+      for j=1:nvars
+        temp_vd(:,i,j) = vd(:,((j-1) + ((i-1)*nvars)+1));
+      end
+    end
+    shocks = temp_vd;
+  end
+  
+  for i=1:nvars
+    shock_names{i} = endo_names(i,:);
+    series_names{i} = endo_names(i,:);
+  end
+
+  if nargin < 2
+    title_ = '';
+  end
+
+  x = [1:T]; plot_dates = 0;
+	
+	data = 0;
+	steady = 0;
+	
+	colors = [ .1 .1 .75
+			 .8 0 0
+			 1 .7 .25
+			 1 1 0	
+			 .5 1 .5 
+			 .7 .7 .1
+			 .5 .6 .2
+			 .1 .5 .1];
+	
+	% overide the defaults with optional inputs
+	for i=1:length(varargin)
+		if strcmpi(varargin{i},'data')
+			data = varargin{i+1};
+		elseif strcmpi(varargin{i},'steady')
+			steady = varargin{i+1};
+		elseif strcmpi(varargin{i},'shock_names')
+			shock_names = varargin{i+1};
+		elseif strcmpi(varargin{i},'series_names')
+			series_names = varargin{i+1};
+		elseif strcmpi(varargin{i}, 'dates')
+			x = varargin{i+1}; plot_dates = 1;
+		elseif strcmpi(varargin{i},'colors')
+			colors = varargin{i+1};
+		end
+	end
+	
+	% add an extra period to the time series
+	x(T+1) = x(T) + (x(T) - x(T-1));
+
+	
+	figure('Name',title_)
+	for k=1:K
+		% Go through each series
+		subplot(K,1,k);
+		sshocks = shocks(:,k,:);
+		hold on
+		% plot the stacked shocks