Monte Carlo Filtering
============================

.. index:: filtering

Description
----------------

Monte Carlo filtering is a method for analyzing the sensitivity of a model's characteristics to its parameters. The
approach involves generating draws of parameters from their distribution (e.g. prior, posterior, support) and
examining the model's behavior for each draw.

By comparing the empirical marginal distributions of the
draws that lead to desired behavior versus those that do not, key parameters that drive the model's behavior
can be identified.

Monte Carlo filtering is useful for locating the most important parameters of a model that contribute to its behavior.

A key reference in Global sensitivity analysis is :cite:`Ratto2008`.

Quick start examples
----------------------

A very simple example
~~~~~~~~~~~~~~~~~~~~~~~~

.. code-block:: matlab
   :linenos:

	x=randn(2^13,8);
	
	y=x(:,1)+x(:,2).^2+x(:,1).*x(:,2);
	
	check_behavior=y<0;
	
	obj=mcf(check_behavior(:).',x.')
	%%
	cdf_plot(obj)
	%%
	scatter(obj)
	%%
	correlation_patterns_plot(obj,[],'behave')

In this small example, we have a function y=y(x1,x2) and we interested in finding the combinations of parameters/covariates
x1, x2, x3, x4, x5, x6, x7 and x8 for which y<0. In this controlled example, clearly, we know a priori that the only
covariates that matter for the result are x1 and x2. But in general we may not know that in advance.

We first sample values for all the covariates and for each combination (x1(i), x2(i), x3(i), x4(i), x5(i), x6(i), x7(i), x8(i)),
for i=1,2,...,2^13, we compute y(i) and set a flag telling us whether y(i) is negative for combination i.

Then we proceed to investigate whether there are some patterns among the different covariates both in situations where y<0
and in situations where y>=0.

A simple DSGE example
~~~~~~~~~~~~~~~~~~~~~~~~

.. code-block:: matlab
   :linenos:

	%% A DSGE model
	
	simpleModel={
	    'model:someModel'
	    '%% model without parameters and shocks'
	    '@endogenous X1, X2, X3'
	    '@exogenous E1, E2, E3'
	    '@parameters a1, a2, a3, vol_tp_1_2, vol_tp_2_1'
	    '@parameters(vol,2) s1, s2, s3'
	    '@model'
	    'X1=a1*X1{-1}+s1*E1;'
	    'X2=a2*X2{-1}+s2*E2;'
	    'X3=a3*X3{-1}+s3*E3;'
	    };
	
	%% The parameters to vary and their ranges
	
	plist={'a1', 'a2', 'a3', ...              
	    'vol_tp_1_2', 'vol_tp_2_1',...        
	    's1_vol_1', 's2_vol_1', 's3_vol_1',...
	    's1_vol_2', 's2_vol_2', 's3_vol_2'    
	    };                                    
		
	bounds=[0.5,0.9;0.5,0.9;0.5,0.9;0,1;0,1;0,0.5;0,0.5;0,0.5;0.5,2;0.5,2;0.5,2];
	
	%% rise the model and set the behavior to investigate in a function handle
	
	m=rise(simpleModel);

	% To do : write a simple function handle check that the model solves or that some correlations match
	% Behavior = @...

	%% run the behavior through the mcf

	obj=mcf(Behavior,nsim_or_draws,bounds(:,1),bounds(:,2),plist,'latin_hypercubes');


	%% smirnov test for equality of distributions
	close all
	% hfig=utils.plot.multiple(@(x)cdf_plot(obj,x),parameter_names,...
	%     'Smirnov test of equality of distributions',3,3)
	% 
	cdf_plot(obj)
	
	pause
	
	%% correlation patterns in the behavior sample
	close all
	
	correlation_patterns_plot(obj,[],'')
	
	pause
	
	%% correlation patterns in the behavior sample
	close all
	
	correlation_patterns_plot(obj,[],'behave')
	
	pause
	
	%% correlation patterns in the behavior sample
	close all
	
	correlation_patterns_plot(obj,[],'non-behave')
	
	pause
	
	%% scatter plot of the significant correlations
	close all
	
	scatter(obj,[],'',1/100)
	
	pause
	%% scatter plot of the significant correlations
	close all
	 
	scatter(obj,[],'non-behave')

This second example transposes the workflow above to a DSGE model with
regime-switching shock volatilities. The model has three AR(1)
endogenous variables ``X1``, ``X2``, ``X3``, three shocks ``E1``,
``E2``, ``E3``, and two volatility regimes (``vol``) governed by the
transition probabilities ``vol_tp_1_2`` and ``vol_tp_2_1``. The
standard deviations ``s1``, ``s2``, ``s3`` are declared *switching*
over the ``vol`` chain, so each one carries two values in the
parameter table (``s1_vol_1`` / ``s1_vol_2``, etc.).

The parameter list ``plist`` collects every parameter the analyst
wants to vary: the AR(1) persistences ``a1``-``a3``, the two
transition probabilities, and both regime values of each volatility.
``bounds`` gives the lower and upper sampling range for each.

The placeholder ``Behavior`` would be a function handle of the form
``@(x) someTest(x)`` that returns ``true`` when the sampled parameter
vector ``x`` produces the model behaviour we want to single out --
for instance, that the model solves, that a particular impulse
response falls in a target range, or that a moment condition is met.
The following section (*Support facility for filtration*) shows
``mcfprop`` as a convenient builder for that handle.

With a ``Behavior`` handle in place, ``mcf`` runs the same battery
as in the previous example: ``cdf_plot`` reports the marginal
distribution of each parameter conditional on behaviour vs.
non-behaviour; ``correlation_patterns_plot`` highlights pairwise
patterns within the behave / non-behave / full samples; ``scatter``
shows the joint distribution of the most informative pairs.

Details on the mcf class
--------------------------

.. toctree::
   :maxdepth: 2
   :caption: Contents:

   monte_carlo_filtration

Support facility for filtration : the mcfprop class
-----------------------------------------------------

.. toctree::
   :maxdepth: 2
   :caption: Contents:

   monte_carlo_filtration_support

An example
~~~~~~~~~~~~

.. code-block:: matlab
   :linenos:
	
	%% construct the support object
	
	m=rise(simpleModel);
	
	self = mcfprop(m, plist);
	
	%% Add a few properties
	self = self.add('solve'); % the model solves
	self = self.add({'corr{X2,X2,3:8,2}',@(x)x>0.01}); % correlation or order 3 in regime 2
	self = self.add({'corr{X3,X3,0,2}',@(x)abs(x-1)<1e-9}); % correlation or order 0 in regime 2
	self = self.add({'irf{X1,E1}',@(x)x>0.01}); % irf in period 1
	self = self.add({'irf{X1,E1,3}',@(x)x>0.01}); % irf in period 3
	self = self.add({'irf{X2,E2,3:5,2}',@(x)all(x>1e-3)}); % irf in period 3 to 5 in regime 2
	self = self.add({'Tz{X1,E1}', @(x) x > 0}); % regime 1
	self = self.add({'Tz{X1,@sig,2}', @(x) abs(x)<1e-8}); % regime 1
	
	%% Create a function handle for the behavior
	
	n=size(bounds,1);
	
	propfun = create_handle(self);
	
	%% Evaluate one parameter vector
	
	x = bounds(:,1)+ rand(n,1).*(bounds(:,2)-bounds(:,1));
	
	eval_ = propfun(x);

``mcfprop`` is the builder side of ``mcf``: it lets the user assemble
a behaviour-checking function from named model properties without
writing it by hand. Each call to ``self.add(...)`` registers one
property to evaluate. The supported property names mirror RISE's
post-solve accessors:

- ``'solve'`` -- the model solved successfully at the candidate
  parameter vector.
- ``'corr{Xi,Xj,k:l,r}'`` -- the autocovariance / cross-covariance
  of order ``k`` (or range ``k:l``) between variables ``Xi`` and
  ``Xj``, in regime ``r``.
- ``'irf{Xi,Ej,k:l,r}'`` -- the impulse response of ``Xi`` to shock
  ``Ej`` at horizon ``k`` (or range), in regime ``r``.
- ``'Tz{Xi,Ej,r}'`` / ``'Tz{Xi,@sig,r}'`` -- direct entries of the
  policy function in regime ``r`` (impact of shock ``Ej``, or of the
  perturbation parameter ``@sig``).

The optional second element of the cell is a predicate ``@(x) bool``
that the value of the property must satisfy for the parameter vector
to count as behaving. Without a predicate the property's value is
simply recorded.

``create_handle(self)`` collapses the entire menu into a single
callable ``propfun(x)`` whose return value is the joint truth of all
predicates. That handle is exactly what the ``Behavior`` slot of
``mcf`` expects -- closing the loop from a higher-level statement of
"model properties of interest" to a low-level boolean signal driving
the sensitivity analysis.

.. hint:: With the handle/behavior in hand we can go ahead and pass it to the mcf class as above