\documentstyle[aps]{revtex}
%\def\btt#1{{\tt$\backslash$#1}}
\begin{document}
\title{Metropolis Path Integral 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 4}
\vskip 2truecm

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

Path Integral Monte Carlo Simulator (PIMCS) is used to generate complete histories of 
correlated instrument movements which are described by a Markov process. It can be
easily extended to treat pseudo-Markov processes.
A detailed description of the underlying theory can be found in the accompanying document.

Let ${\bf X}_t$ denote a vector of observable rate values (or instruments)
 at time $t$ which has $N$ components. 
It is assumed that
all $N$ rates are correlated. 
The Markov process is completely defined by the conditional probability of reaching
vector ${\bf X}_{t+\Delta t}$ at time $t+\Delta t$ given a vector ${\bf X}_t$ at time $t$:
\begin{equation}
P({\bf X}_{t+\Delta t}, t+\Delta t \vert {\bf X}_t , t)
\label{prob}
\end{equation}
In the limit of infinitesimal $\Delta t$, the conditional probability becomes the Markov
propagator.

For example, the stochastic differential equation which is the basis of derivative pricing
formulas:
\begin{equation}
d{\bf X} (t) = {\bf\alpha} ({\bf X},t)dt + {\bf \sigma}^2 ({\bf X},t){\bf Z}dt
\label{ito1}
\end{equation}
where
${\bf Z}$ is a random vector drawn from $N(0,1)$ can be modeled as a Markov process with the
propagator:
\begin{equation}
P({\bf X}_{t+\Delta t}, t+\Delta t \vert {\bf X}_t , t) \propto e^{-{1\over2} ({  {{\bf X}_{t+\Delta t} - {\bf X}_t - {\bf\alpha} ({\bf X}_t , t) \Delta t} \over
   { {\bf \sigma}^2 ({\bf X}_t , t) \Delta t} })}
\label{markov}
\end{equation}

PIMCS uses conditional probability function and Metropolis algorithm to generate
rate histories. Probability distribution of a rate history is just a product of conditional
probabilities for all time slice of the history.
There are many ways one can generate histories given a conditional probability, but the unique feature
of the Metropolis method is that it can be applied to any probability function, 
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 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 of histories. 
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 
configurations 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, path or configuration will be used interchangeably
to denote rate histories which are generated in the course of a Monte Carlo simulation.
In other words, each of these terms denotes a snapshot of values of all rates for all time steps.


\section{The Path Integral 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 PIMCS are requested through a 4 step process:

\begin{enumerate}
\item {\bf Creation of the data structure which will contain the paths}. 
      A call to
      the routine {\bf CreatePath()} will create 
      a {\bf PATH}  data structure which will be used to store values of all rates for all time steps.
      The pointer to this data structure is
      returned to the caller:
\begin{verbatim}

PATH *Rates;        /* Pointer to the path container               */
int NumTimeSteps;   /* Number of of time steps in a path.          */ 
int NumRates;       /* Number of rates or instruments, 
                       i.e. dimensionality of the rate vector      */
int ProcessOrder;   /* Size of the extra buffer to store initial 
                       lags of rate histories. This buffer is not 
                       used by the PIMCS for true Markov processes, 
                       but may be used for pseudo-Markov processes */

Rates  = CreatePath
                  (
                  NumTimeSteps,
                  NumRates,
                  ProcessOrder
                  )
\end{verbatim}
\item {\bf Creation of the data structure which contains the description of the conditional
probability }. Parameters which define the required probability 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 *ProcessParams;          /* Pointer to the PROCESS data structure, which 
                                     contains parameters of the conditional 
                                     probability                                  */
int NumTimeSteps;
int NumRates;
char *ProcessDescriptionFile;    /* Name of the input file which contains
                                    parameters of the conditional probability     */

ProcessParams   =   SetProcessParams
                                   (
                                   NumTimeSteps,
                                   NumRates,
                                   ProcessDescriptionFile
                                   )
\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
paths: 
\begin{verbatim}

PATH *Rates
PROCESS *ProcessParameters;

GetPathSample
         (
         Rates,
         ProcessParams
         );
\end{verbatim}
\end{enumerate}

PIMCS starts from an initial path , generates a Markov chain of 
paths , 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 path reached in the Metropolis Markov chain
is returned as the proxy rate history}. If requested, analysis can be done as the paths are
being generated. For example, in the test code, second moments of the
rates 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 Path Integral Monte Carlo Kernel}
\label{sec:kernel}

The {\bf GetPathSample()} routine is a wrapper around the Metropolis Monte Carlo kernel,
\begin{verbatim}
static int MetropolisMonteCarlo
                              (
                              PATH *Rates,
                              PROCESS *ProcessParams
                              )

{
  int mb;
  int mcs;
  int  Naccepted=0;
  int NumRates = HowManyRates(Rates);
  int NumMeasBlocks   = HowManyMeasBlocks();
  int NumRelaxBlocks  = HowManyRelaxBlocks();
  int MeasBlockSteps  = HowBigMeasBlock();

  SetAllPathsToZero(Rates);

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

  while(!ImDone()) {
    for(mb=1;mb<=NumMeasBlocks; mb++) {
      Naccepted += MetropolisUpdate
                                  (
                                  MeasBlockSteps,
                                  Rates,
                                  ProcessParams
                                  );
    }
    CycleDone();
  }

  return Naccepted;

} /* end of MonteCarlo() driver */

   
\end{verbatim}

The structure of the Path Integral Metropolis driver is identical to the other Metropolis
driver, and is described in the documentation for the Metropolis Monte Carlo Simulator
(technical note \# 3).

There is an important difference in the functionality of the {\bf MetropolisUpdate()} routine,
although the interface is identical in both drivers. The fundamental object in PIMCS is a path.
{\bf MetropolisUpdate()} generates a block of successive paths along the Metropolis Markov chain.
In order to use uncorrelated configurations, only the last path in a block is 
recorded in the {\bf PATH *Rates} container, and used for analytics. 

Given a current path, the trial path is obtained by adding a perturbation to a single rate on a
single time slice. Due to the Markov property, probability of the trial path depends on the
rate values at the chosen time slice, as well as the previous time slice and the next time
slice. Metropolis test is performed using the ratio of probabilities of the trial path and the
current path. If the trial path is accepted, it becomes the current path. This procedure is
iterated for all rates and all time slices sequentially. When an update is attempted for
all rates on all time slices, one Monte Carlo step is completed. If there are {\bf MeasBlockSteps}
in a block, {\bf MetropolisUpdate()} will perform {\bf MeasBlockSteps} Monte Carlo steps.
The final path in a block is stored in {\bf Rates} and the control  is returned to the driver.
The driver will iterate {\bf MetropolisUpdate()}  {\bf NumRelaxBlocks} times for relaxation and
then {\bf NumMeasBlocks} times for analytics calculations. If only proxy rate generation is
required then {\bf NumMeasBlocks} may be set to $1$. 


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

{\bf PROCESS} data structure contains fields which identify the conditional probability function
as well as its parameters. {\bf MetropolisUpdate()} will obtain the address of the conditional probability
function and its parameters by issuing appropriate inquiries to the module which manages {\bf PROCESS}
data structure.

Conditional probability functions may be  accessed using two interfaces:
\begin{enumerate}
\item Interface to the conditional probability function which {\bf bundles 
      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 path updates for all rates on a single time slice
      to the conditional probability function. The probability function  computes probabilities and
      performs Metropolis acceptance/rejection test. If a proposed update is rejected, the
      corresponding element of the {\bf TrialPerturbations} array is set to $0$. 
      This interface allows more efficient untilization of vector units, if available.
\begin{verbatim}
int NumRates;
float *TrialPerturbations;  /* array of proposed updates for all rates on 
                               the present time slice                    */
float *RatesPres;           /* rate values on the present time slice     */
float *RatesPrev;           /* rate values on the previous time slice    */
float *RatesNext;           /* rate values on the next time slice        */
float **Couplings;          /* Parameters of the conditional probability 
                               function. 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; the interpretation
                               of the matrix is up to the 
                               ConditionalProbabilityFunction()          */ 
                    
int ConditionalProbabilityFunction
                                 (
                                 NumRates,
                                 TrialPerturbations,
                                 RatesPres,
                                 RatesPrev,
                                 RatesNext,
                                 Couplings
                                 )
\end{verbatim}
Probability function
returns the number of accepted Metropolis trials. It also has a side effect of changing the
{\bf TrialPerturbations} array. If a proposed update is rejected, the 
corresponding element of the {\bf TrialPerturbations} array is set to $0$, otherwise it is left as is.
When the
control is returned from the {\bf ConditionalProbabilityFunction()} to {\bf MetropolisUpdate()},
{\bf TrialPerturbations} are added to the rate values on the present time slice ({\bf RatesPres}). 

\item Interface to the {\bf probability function} which {\bf only performs transition probability 
      computation}.  
      This interface passes a proposed update of a single rate 
      to the conditional probability function, which computes the probability of acceptance of the trial
      path.
      Metropolis acceptance/rejection test and updates of the current path are performed
      within MetropolisUpdate(). This interface provides more flexibility within 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 *RatesPres;
float *RatesPrev;
float *RatesNext;
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 ConditionalProbabilityFunction
                                   (
                                   NumRates,
                                   TrialIndex,
                                   shift,
                                   RatesPres,
                                   RatesPrev,
                                   RatesNext,
                                   Couplings
                                   )

\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)
7
\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 paths from the path 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 on every time step.
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
path, 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
paths 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 conditional probability function (1= Gaussian, 3=Student(n=3), ...)
1
Standard deviations of rates (used to set the scale of Metropolis updates)
1 (float)
2 (float)
.
.
.
m (float)
Parameters of the conditional probability function (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 conditional probability function. 
After that, standard deviations of each of the component of
the rate vector are listed, one per line. Each line gives the rate 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 rates is irrelevant.
Next,
the parameters of the conditional probability function are given as a matrix. 
Number of rows is equal to the
number of rates ($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 time steps ina a history (temporal extent of a path)
\item  Dimensionality of the rate vector  (number of rates )
\end{enumerate}
and addresses of the following data structures, which must be prepared before the simulator
is called:
\begin{enumerate}
\item {\bf PATH} *RateHistories
\item {\bf PROCESS}  *ConditionalProbabilityDescriptor
\end{enumerate}

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

The Monte Carlo Simulator will generate paths of specified dimensionality (number of rates, number of
time steps)
and place them in {\bf RateHistories}.

\end{document}
