\documentstyle[aps]{revtex}
%\def\btt#1{{\tt$\backslash$#1}}
\begin{document}
\title{Metropolis Monte Carlo Simulator\\
}
\vskip 1truecm
\renewcommand{\thefootnote}{*}
\author{Miloje S. Makivi\'c\footnote{
Northeast Parallel Architectures Center (NPAC), Syracuse University,\\
111 College Place, Syracuse, NY 13244-4100\\ 
{\bf e-mail}: miloje@npac.syr.edu\\ 
{\bf Mosaic} http://www.npac.syr.edu/users/miloje/}
}
\date{\today}
\maketitle
\vskip 2truecm
\centerline{NPAC/Sailfish  Technical Note 3}
\vskip 2truecm

\section{Overview}
\label{sec:Overview}

Monte Carlo Simulator is used to generate vectors which are sampled from a specified probability
distribution. As an example of an application of the simulator, consider the
module which generates vector auto-regressive time series of order $k$:
\begin{equation}
{\bf X} (t) = {\bf\alpha} + {\bf A_1}{\bf X}(t-1) + {\bf A_2} {\bf X} (t-2) +\dots + {\bf A_k} {\bf X} (t-k) + {\bf U} (t)
\label{var}
\end{equation}
with a vector auto-regressive conditional heteroskedastic noise of order $m$:
\begin{equation}
{\bf U} (t) = ({\bf\xi} + {\bf B_1} {\bf U}^2 (t-1) + {\bf B_2} {\bf U}^2 (t-2) +\dots + {\bf B_m} {\bf U}^2 (t-m))^{1/2} \epsilon(t)
\label{arch}
\end{equation}
A detailed description of the implementation of this module is given in the accompanying document.
${\bf X} (t)$ is vector of observed rates at time $t$ which has $N$ components. It is assumed that
all $N$ rates are correlated. Rate levels are given by ${\bf\alpha}$.
${\bf X} (t)$ is regressed on $k$ previous time steps via $k$ matrices
${\bf A_p},p=1\dots k$ which are assumed to be time independent. ${\bf U} (t)$ is an $N$-component
noise vector which is also regressed onto $m$ previous time steps via ${\bf B_p},p=1\dots m$ matrices.
Noise levels are given by vector ${\bf\xi}$.

Noise vector contains $N$-component {\bf disturbance vector} $\epsilon(t)$ which {\bf is sampled} from an
arbitrary distribution {\bf using the Metropolis Monte Carlo Simulator}. 
This would be a typical Sailfish application of the Monte Carlo Simulator, but it is certainly not limited
to time series modeling. 

There are many ways one can generate vector samples from specific distributions, but the unique feature
of the Metropolis method is that it can be applied to any distribution and it is going to work efficiently
on multi-dimensional problems. Given that the dimensionality of Sailfish applications is of the order
of $10^2$ (i.e. at least $100$ coupled instruments) and required flexibility to quickly implement
different noise processes, Metropolis seems to be the method of choice.

Metropolis works by constructing a Markov chain in
the phase space of the problem. This is a great advantage of the method, since one needs only ratios of
probabilities of phase space points to construct the transition probability between two successive
states along the Markov chain. The invariant distribution of the Markov chain is the desired probability
distribution. The generality of the method comes with a price: successive states generated in the course
of the simulation
are correlated, which is not the case with other Monte Carlo methods. To obtain uncorrelated samples one
usually retains only a fraction of generated configurations which are sufficiently far apart along the
Markov chain. This implies that Metropolis method is usually slower than other methods, but I believe
that its flexibility offers far bigger advantages. Another reason is that eventually Sailfish applications
will use Monte Carlo Simulator to perform analysis on generated instrument histories. Depending on the
nature of the analysis, its computational cost will most likely be large enough to amortize the
cost of Metropolis or any other Monte Carlo algorithm.

In the rest of the text, vector sample or configuration will be used interchangeably
to denote random vectors which are generated in the course of a Monte Carlo simulation.


\section{The Monte Carlo Simulator Interface}
\label{sec:interface}

This document
describes only the mechanics of the simulator. For a detailed
description of relevant data structures, please see the accompanying document.
Services of the Monte Carlo Simulator are requested through a 4 step process:

\begin{enumerate}
\item {\bf Creation of the data structure which will contain the vector samples}. For example, the
      ARCH noise module requires a disturbance vector for each time step of the series. A call to
      the routine {\bf CreatePath()} will create 
      a {\bf PATH}  data structure which will be used to store disturbance vectors for all time steps.
      The pointer to this data structure is
      returned to the caller:
\begin{verbatim}

PATH *Disturbance;  /* Pointer to the container for vector samples              */
int NumTimeSteps;   /* Number of vector samples. In the context of time series
                       analysis, this is the number of time steps in the series */
int NumRates;       /* Dimensionality of the vector sample. In the context
                       of time series analysis, this is the number of rates or
                       instruments                                              */
int ProcessOrder;   /* The order of the time series. This variable determines 
                       the size of the extra buffer to store initial lags of the
                       time series. This buffer is not used by the Monte Carlo 
                       Simulator, but it is used by other modules               */

Disturbance = CreatePath
                       (
                       NumTimeSteps,
                       NumRates,
                       ProcessOrder
                       )
\end{verbatim}
\item {\bf Creation of the data structure which contains the description of the sample 
probability distribution}. Parameters which define the required distribution are stored
in a flat ASCII file. Routine {\bf SetProcessParams()} will open the file, read the parameters,
create a {\bf PROCESS} data structure where the parameters are stored and return a pointer to the data
structure to the caller:
\begin{verbatim}

PROCESS *DisturbanceParams;       /* Pointer to the PROCESS data structure, which 
                                     contains the parameters of the probability
                                     distribution of vector samples               */
int NumTimeSteps;
int NumRates;
char *DisturbanceDescriptionFile; /* Name of the input file which contains
                                     the parameters of the probability
                                     distribution of vector samples               */

DisturbanceParams = SetProcessParams
                                   (
                                   NumTimeSteps,
                                   NumRates,
                                   DisturbanceDescriptionFile
                                   )
\end{verbatim}
\item {\bf Setup of the control structure for Metropolis simulation}. Routine
{\bf SetupMetropolisMCSimulation()}  will open an input file which contains parameters related to
the mechanics of the Metropolis method, read the parameters and configure the Metropolis driver
accordingly. For example, the input file will specify number of relaxation steps, number of steps
between recorded configurations (which is used to eliminate correlations between configurations),
initial seed for the random number generator, total length of the simulation, etc. Metropolis
algorithm specific parameters are stored in a MONTECARLO data structure, which is not visible
to the rest of the application.
\begin{verbatim}
char *MC_description_file; /* input file Metropolis algorithm specific parameters */
SetupMetropolisMCSimulation(MC_description_file)
\end{verbatim}
\item {\bf Issuance of a request to do the simulation} and fill the supplied data structure with
the vector samples. For example, ARCH noise module will obtain disturbance vectors by calling:
\begin{verbatim}

PATH *Disturbance;
PROCESS *DisturbanceParameters;

GetDisturbanceSamples
         (
         Disturbance,
         DisturbanceParams
         );
\end{verbatim}
\end{enumerate}

Monte Carlo Simulator starts from a zero vector, generates a Markov chain of 
vector samples , which is, after the specified number of relaxation steps, drawn
from the requested distribution. For the purpose of time series analysis and proxy rate
histories generation, the {\bf final configuration reached in the Markov chain
is returned as the disturbance vector sample}. If requested, analysis can be done using
generated disturbance configurations. For example, in the test code, second moments of the
disturbance are computed during the simulation to verify its accuracy. Further analytics
modules (for example, option pricing) will be added later on as a part of the {\bf Analysis
Module}.


\section{The Metropolis Monte Carlo Kernel}
\label{sec:kernel}

The {\bf GetDisturbanceSamples()} routine is a wrapper around the Metropolis Monte Carlo kernel,
\begin{verbatim}
static int MetropolisMonteCarlo
                              (
                              PATH *Disturbance,
                              PROCESS *DisturbanceParams
                              )
   
{
  int NumMeasBlocks;
  int MeasBlockSteps;
  int NumRelaxBlocks;
  int mb;
  int mcs;
  int Naccepted=0;
  int NumRates = HowManyRates(Disturbance);

  NumMeasBlocks   = HowManyMeasBlocks();
  NumRelaxBlocks  = HowManyRelaxBlocks();
  MeasBlockSteps  = HowBigMeasBlock();

  SetAllPathsToZero(Disturbance);  /* begin from a zero vector                                */

  /* Samples are asymptotically drawn from the desired probability distribution. In order 
     to eliminate the influence of the initial configuration, a number of samples 
     (= NumRelaxBlocks) at the beginning of the simulation is discarded                       */ 

  SetRelaxFlag(NOT_RELAXED);
  if(!IsRelaxed()) for(mb=1;mb<=NumRelaxBlocks; mb++)
            Naccepted += MetropolisUpdate(MeasBlockSteps, Disturbance, DisturbanceParams);
  SetRelaxFlag(RELAXED);

  /* After the relaxation is completed, samples are drawn from the specified probability
     distribution. Successive configurations generated in the course of the simulation
     are correlated. Because of this, configurations are recorded with a user specified
     frequency (=MeasBlockSize). Simulation is broken into blocks which contain MeasBlockSize
     successive configurations. Configurations within a block are generated using the
     MetropolisUpdate() routine.  Only the last configuration in a block is recorded in 
     the Disturbance structure and used for analysis. MeasBlockSize is chosen to be larger 
     than the correlation length along the Metropolis Markov chain. Total number of 
     configuration blocks is equal to the number of calls to MetropolisUpdate(), which is
     in turn equal to (NumMeasBlocks + NumRelaxBlocks).                                       */
     

  while(!ImDone()) {

    for(mb=1;mb<=NumMeasBlocks; mb++) {
      Naccepted += MetropolisUpdate(MeasBlockSteps,  Disturbance, DisturbanceParams);
    }
    CycleDone();     /* Simulation can be further structured into a number of cycles.
                        After a cycle is completed, the complete state of the Monte 
                        Carlo Simulater can be dumped onto a disk for checkpointing/
                        restarting. This is a good strategy to prevent data loss in 
                        the event of system failure for long batch jobs. At the end
                        of a cycle, analytics data can be also dumped into files for
                        later processing. For proxy rate generation, number of cycles
                        is always 1.                                                         */

  }

  return Naccepted;

} /* end of MonteCarlo kernel */
\end{verbatim}

In a Metropolis algorithm, one generates a sequence of configurations or samples, which are
asymptotically drawn from the specified probability distribution. Given a present configuration,
one first generates a trial new configuration. The trial configuration is accepted with probability
which is equal to the ratio of probabilities of the trial configuration and the current 
configuration. If the trial configuration
is accepted it becomes the current configuration and the procedure is iterated as many times as desired.
The computation of probabilities and acceptance/rejection test are done within
{\bf MetropolisUpdate()}. {\bf MetropolisUpdate()} requires that a {\bf ProbabilityFunction()} be
identified in the {\bf PROCESS} data structure.




\section{Probability Functions}
\label{sec:prob}

{\bf PROCESS} data structure contains fields which identify the sample probability function
as well as its parameters. {\bf MetropolisUpdate()} will obtain the address of the probability
function and its parameters by issuing appropriate inquiries to the module which manages {\bf PROCESS}
data structure.
Probability functions may be  accessed using two interfaces:
\begin{enumerate}
\item Interface to the probability function which {\bf bundles 
      transition probability computation and Metropolis test}. This interface is chosen by
      defining the BUNDLE\_AR\_WITH\_PROB preprocessor switch. This interface passes a whole
      vector of proposed updates to the probability function, which computes probabilities,
      performs Metropolis acceptance/rejection test and updates the current configuration.
      This interface allows more efficient untilization of vector units, if available.
\begin{verbatim}

int NumRates;
float  *TrialPerturbationsAroundCurrentVectorSample;
        /* In a Metropolis algorithm, one generates a trial new configuration,
           given the present one. This array contains proposed updates for all
           components of the current sample (i.e. its dimensionality is
           NumRates). In the TransitionProbability() routine, each component
           of the current sample will be added the respective perturbation value,
           and the new trial value will be accepted or rejected using the Metropolis   
           test. If the trial value is accepted, it will overwrite the old value
           in the current configuration array                                      */
float  *CurrentVectorSample;
        /* Array which stores current sample (configuration)                       */
float  **ParametersOfTheSampleProbabilityFunction;
        /* Many probability functions can be parametrized by a matrix, which is why
           it is chosen as the format for parameter arguments. However, one may
           pass any parameter structure using this format, and the interpretation
           of the matrix is entirely up to the TransitionProbability() routine     */ 
           

int TransitionProbability 
                        (
                        NumRates , 
                        TrialPerturbationsAroundCurrentVectorSample, 
                        CurrentVectorSample,
                        ParametersOfTheSampleProbabilityFunction
                        );
\end{verbatim}
Probability function
returns the number of accepted Metropolis trials. It also has a side effect of updating the
{\bf CurrentVectorSample}.  After a call to {\bf TransitionProbability()}, {\bf CurrentVectorSample}
contains a new configuration of the Metropolis Markov chain. 

\item Interface to the {\bf probability function} which {\bf only performs transition probability 
      computation}.  
      This interface passes a proposed update of a single component of the sample vector
      to the probability function, which computes the probability of acceptance of the trial
      configuration.
      Metropolis acceptance/rejection test and updates of the current configuration are performed
      within MetropolisUpdate(). This interface provides more flexibility with the updating procedure,
      which may be useful if symmetries of the transition probability can be used to improve
      sampling efficiency. Also, for the exponential family of probabilities it may reduce the
      number of parameters that are passed on the stack. 
\begin{verbatim}

int NumRates;
int TrialIndex;  /* the index of the vector component which is updated */
float shift;     /* the value of the proposed update                   */
float *CurrentVectorSample; 
float **Coupling /* the parameters of the probability function  
                    Depending on the probability function, if the update
                    of a single vector component does not require all
                    parameters, then this parameter may be reduced to
                    a single pointer. This is, for example, the case
                    with Gaussian distribution.                       */
float TransitionProbability
                          (
                          NumRates,
                          TrialIndex,  
                          shift,
                          CurrentVectorSample, 
                          Coupling 
                          )
\end{verbatim}
Probability function returns the probability of the proposed update. It is the responsibility of
the caller (i.e. {\bf MetropolisUpdate()} to accept/reject and perform the update.

\end{enumerate}

\section{Input}
\label{sec:input}
Monte Carlo Simulator requires two input ASCII files:
\begin{enumerate}
\item Metropolis Simulation Description File 
\item Probability Distribution Description File
\end{enumerate}
Metropolis Simulation Description File has the following format:
\begin{verbatim}
Number of Metropolis cycles ( should be 1, unless check-pointing is done)
1
Initial seed for the random number generator
4567
Desired number of uncorrelated MC configurations (= NumMeasBlocks)
500
Estimated number of MC steps between 2 uncorrelated configurations (= MeasBlockSize)
10
Number of uncorrelated steps to be used for relaxation(= NumRelaxBlocks)
10
\end{verbatim}
Each parameter is on a separate line. Each parameter line is preceeded by a comment line.
In this example,
we estimate that correlation length along the Markov chain is 10 steps or less.
We want to generate 500 independent samples from the probability distribution.
This requires generating 500*10 Monte Carlo steps. A Monte Carlo step is defined
as a sweep during which one tries to update every component of the rate vector.
Also, to assure that we are indeed drawing from the desired probability distribution
and not from some transient distribution which still has the memory of the initial
configuration, we run the simulation until we generate  7 uncorrelated configurations
(i.e. 7*10 Monte Carlo steps). After that, we may start using the newly generated
configurations for computation. As it stands now, the Monte Carlo driver does not
do check-pointing and restarting from a file dump, but all the facilities for that are
provided. This is useful to prevent long simulations to be completely wasted by
a system crash. Number of MC cycles is equal to the number of checkpoints and
restarts. For the present purpose of generating proxy rate histories
it should be set to 1, but it is not a mistake
to set it to any other number.

Probability Distribution Description File has the following format:
\begin{verbatim}
Index of the probability distribution function (1= Gaussian, 3=Student(n=3), ...)
1
Standard deviations of disturbance vector (used to set the scale of Metropolis updates)
1 (float)
2 (float)
.
.
.
m (float)
Parameters of the disturbance probability distribution (see below for description of format)
1 N1
1 1 (float)
1 2 (float)
.
.
1 N1 (float)
2 N2
2 1 (float)
2 2 (float)
.
.
.
2 N2 (float)
.
.
.
m Nm 
m 1 (float)
m 2 (float)
.
.
.
m Nm (float)
\end{verbatim}

Parameters are separated by comment lines. The first parameter is an index which specifies
the requested probability function. After that, standard deviations of each of the component of
the sample vector are listed, one per line. Each line gives the component index (from $1$ to $m$)
and its standard deviation (float). Only a 
rough estimate is needed since it is used to establish a scale on which trial configurations
in the Metropolis algorithm differ from the current one.
Since the indices are explicitly listed, the ordering of components is irrelevant.
Next,
the parameters of the distribution function are given as a matrix. Number of rows is equal to the
number of vector components ($m$). Number of columns may
be variable, which can be used to reduce storage requirements when probability function has
symmetries. Because of that, parameters from each row are preceded by a line which says what is
the row number (k) and the number of parameters in that row ($N_k$). Along with each parameter, its
row and column is passed. Because of this, the order in which parameters are listed in the file
is irrelevant. 
The interpretation of the input matrix is entirely up to the corresponding  probability function.
In the case of Gaussian and Students distribution which have quadratic forms $\sum A_{ij} X_i X_j$,
these are rows of the {\bf A} matrix.

In addition to the two input files, Monte Carlo Simulator requires following integer
input parameters:
\begin{enumerate}
\item  Number of sample vectors
\item  Dimensionality of sample vectors
\end{enumerate}
and addresses of the following data structures, which must be prepared before the simulator
is called:
\begin{enumerate}
\item {\bf PATH} *SampleContainer
\item {\bf PROCESS}  *SampleProbabilityDescriptor
\end{enumerate}

\section{Output}
\label{sec:output}

The Monte Carlo Simulator will generate the requested number of sample vectors of requested
dimensionality and place them in {\bf SampleContainer}.

\end{document}
