\documentstyle[art10,npac_proto,palatino]{article}

\begin{document}
%\makeNPACtitle\newpage
\begin{project}
\begin{titlepage}
	\projtitle{Risk}
	\doctitle{NPAC Risk Management Package}
		\author{Miloje S. Makivic}
		\date{January 23, 1995}
		\version{0.1}
		\summary{Draft}
	\id{risk}
	\docauthor{M~S~Makivic}
	\doccontributor{H~W~Yau}
	\security{Commercial in Confidence}
\maketitle
	\approval{Geoffrey C. Fox}
\vfill
%	\begin{history}
%		\release{
%			\author{Miloje S. Makivic}
%			\date{August 3, 1995}
%			\version{0.1}
%			\summary{Initial draft}
%			}
%                \release{
%                        \author{Miloje S. Makivic}
%                        \date{November 1, 1995}
%                        \version{0.2}
%                        \summary{Completion of analytics modules}
%                        }
%
%	\end{history}
\end{titlepage}

\pagenumbering{roman}\setcounter{page}{1}
\begin{abstract}
This document reviews major components of the risk management analytical software developed at
Northeast Parallel Architectures Center. Our risk management package is a collection
of flexible and robust modules written in "C" and designed to handle large scale computational 
problems encountered in typical
risk management applications. The present document describes the algorithms and
functionality of the two major modules in the risk management 
package: Monte Carlo Simulator and Estimation Module. 
\end{abstract}

{\parskip=0pt   % hack to patch format of \tableofcontents
\tableofcontents
%\newpage
\medskip\hspace{2.0cm}\hrulefill\hspace{2.0cm}\bigskip
}

\pagenumbering{arabic}\setcounter{page}{1}

\section{Introduction}
Risk management encompasses some of the most compute-intensive tasks in 
financial engineering. The most sophisticated risk management techniques when
applied to large portfolios are bona fide supercomputer applications. We have
developed an array of techniques to tackle the most complex risk management
problems and we are continually expanding the arsenal of analytical tools and
integration techniques. This document is not concerned with integration since
it is application specific, however network centric computing and client server
architecture are featured prominently among our design principles for any application.
The document focuses on Monte Carlo simulations and optimization techniques for
parameter estimation of the
dynamics of risk factors, which are at the core of any risk management application.

\section{Monte Carlo Simulators}

Path Integral Monte Carlo Simulator (PIMCS) is used to generate complete histories of 
correlated risk factor movements which are described by an arbitrary Markov process. 
It can be
easily extended to treat pseudo-Markov and non-Markov processes.
A detailed description of the underlying theory can be found in \cite{art:mm:1995}.

Let ${\bf X}_t$ denote a vector of observable risk factor values (or instruments)
 at time $t$ which has $N$ components. 
It is assumed that
all $N$ factors are correlated. 
The Markov process is completely defined by the conditional probability of reaching
vector ${\bf X}_{t+dt}$ at time $t+dt$ given a vector ${\bf X}_t$ at time $t$:
\begin{equation}
P({\bf X}_{t+dt}, t+dt \vert {\bf X}_t , t)
\label{prob}
\end{equation}
In the limit of infinitesimal $dt$, 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+dt}, t+dt \vert {\bf X}_t , t) \propto e^{-{1\over2} ({  {{\bf X}_{t+dt} - {\bf X}_t - {\bf\alpha} ({\bf X}_t , t) dt} \over
   { {\bf \sigma}^2 ({\bf X}_t , t) dt} })}
\label{markov}
\end{equation}

PIMCS uses conditional probability function and Metropolis algorithm to generate
factor histories. Probability distribution of a factor 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 large scale risk management applications
can be of the order
of thousands and desired flexibility to quickly implement
additional 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. Because of this, Metropolis is not the method of choice for small scale applications.
In the rest of the text, path or configuration will be used interchangeably
to denote risk factor 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 factors for all time steps.

In another typical situation, dynamics of risk factors is described by a time series model
and the purpose of a Monte Carlo Simulator is to generate the noise term.
We have developed a 
Monte Carlo Simulator specifically to generate vectors sampled from an arbitrary
user supplied probability
distribution. The distribution can even be supplied as an experimentally obtained
histogram. 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}
${\bf X} (t)$ is vector of observed factors at time $t$ which has $N$ components. It is assumed that
all $N$ factors are correlated. Factor 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 application of the Monte Carlo Simulator, but it is certainly not limited
to time series modeling. 

%=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
\section{The Laplace Global Optimizer}

Calibration and estimation of parameters governing the dynamics of risk factors is the
centerpiece of any risk management application.
One of the principal methods of time series analysis is the Maximum Likelihood Estimation
of process parameters \cite{book:hl:1991}. 
Model of stochastic factor movement involves either a probability distribution
for noise or an underlying Markov process. To estimate the parameters of this probability
distribution (we will refer to this task as calibration), 
one first constructs the log-likelihood function using an observed sample
of process history and the postulated process probability distribution. 
The set of parameters of the probability distribution which maximizes the log-likelihood
function is the maximum likelihood parameter estimate.
Therefore, the central task in maximum likelihood estimation is the maximization of log-likelihood
function.
The topology of the log-likelihood function can be quite complex, so robust methods for large scale
global optimization must be used. Global optimization methods are required, since there may be
multiple local extrema, depending on the function itself and process sample size.
NPAC developed two global stochastic optimization packages: one is based on Laplace formula and is
similar in spirit to simulated annealing, while the other one uses a stochastic differential
equation solver. 

	The algorithm \cite{art:mp:1970}  finds the global
maximum of a function $f(\vec{x})$ of
$\vec{x}=(x_1,x_2,...,x_i ,...,x_n)$ degrees of freedom. 
This function will be referred to as the objective function.
Location of the maximum is given by the Laplace formula:
\begin{equation}
\vec{x}^* = \lim_{\lambda\to\infty} 
               {
               {\int_D \vec{x}\exp(\lambda f(\vec{x}))d\vec{x}} 
                       \over 
               {\int_D \exp(\lambda f(\vec{x}))d\vec{x}}
               } = \lim_{\lambda\to\infty} <\vec{x}>
\label{laplace}
\end{equation}
The estimation of the first moment $<\vec{x}>$ for each multiplier $\lambda$
is done using the Metropolis Monte Carlo
algorithm. Simulation begins with multiplier $\lambda = 1$ and proceeds
for larger values of multipliers as specified by the user. 
The estimate $<\vec{x}>$ from previous
multiplier is used as the initial point for the simulation with next
multiplier. 

By increasing $\lambda$, vectors $\vec{x}$ generated
in the Metropolis algorithm are more tightly bound to the location of maximum.
Their spread around the average value is a measure of the uncertainty
of the location of maximum. This uncertainty is reduced to $0$ when $\lambda\to\infty$.
Since this limit cannot be achieved in practice, 
user has to specify the desired accuracy of the maximum location.
Two convergence criteria are available:
\begin{itemize}
\item strong convergence
\item weak convergence
\end{itemize}

If the user requests strong convergence criterion,
simulation terminates when the Monte Carlo estimate of
variance for each degree of freedom, $\sigma_i = (<x_i^2 > - <x_i>^2)^{1/2}$ becomes smaller than user
specified accuracy. The accuracy is specified as a small parameter $\epsilon_i$ for each 
degree of freedom $x_i$ in the input
file. The criterion has different forms for degrees of freedom above and
below a user specified threshold $x_i^T$:
\begin{itemize}
\item Strong relative convergence: $\sigma_i \le \epsilon_i \cdot <x_i>$ if
$<x_i> \le x_i^T$
\item Strong absolute convergence: $\sigma_i \le \epsilon_i$ if $<x_i> > x_i^T$
\end{itemize}
These different criteria ensure that time is not wasted trying to converge for near zero parameters,
which are likely to be irrelevant anyway.
The user should ensure that the largest multiplier value $\lambda_{max}$
used in the optimization process is sufficient to obtain the desired accuracy levels.
Otherwise, simulation stops at $\lambda_{max}$ and its accuracy is determined by $\lambda_{max}$.

The weak convergence option for accuracy control is useful for functions with
either high degree of symmetry near the extremum or flat extrema.
If the difference between $\vec{x}^* = <\vec{x}>$ estimates for two 
successive multipliers $\lambda_1 < \lambda_2$ is smaller than the user specified threshold,
the algorithm is considered to have weakly converged. More precisely, user specifies desired
accuracy as an input parameter $\epsilon$.
If $| <x_i (\lambda_2 ) > - <x_i (\lambda_1 )> | \le \epsilon \sigma_i ((\lambda_2 )$, for each
degree of freedom $x_i$, 
weak convergence is achieved.
 It is user's responsibility to
ensure that $\lambda_2$ is sufficiently bigger than $\lambda_1$, otherwise
this convergence criterion is meaningless. In general, this criterion should be
used only if user has some knowledge about function behavior near extremum. 

An important parameter for the Metropolis algorithm is the typical size of a
perturbation (trial update) around the current configuration. For the initial
simulation at $\lambda=1$, this update range is specified by the user. For
all other values of $\lambda$, it is determined using the variance $\sigma_i$ as computed
for the previous value of $\lambda$. This ensures that ranges of trial updates for each
degree of freedom reflect the
scale of variation of the objective function near its maximum. 

In the context of time series analysis, the Laplace package will be used typically
for maximum likelihood estimation, although it can be used for any kind of large
scale optimization problem. 
User has an option to mask the
objective function degrees of freedom which are not variational. This mask is processed
by the Metropolis update routine which ensures that trial moves
are not attempted along the masked coordinates. 

For purposes of time series estimation,
Laplace package computes parameter error estimates using standard procedures.
Log-likelihood function is regarded as a probability function with
Gaussian behavior near the maximum. Parameter error bounds are then obtained
as diagonal elements of the covariance matrix, $(<x_i^2 > - <x_i>^2)^{1/2}$.
Monte Carlo can easily be extended to compute full error covariance matrix,
$(<x_i x_j > - <x_i><x_j >)^{1/2}$, but that would incur quite high
CPU and memory costs.

%=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
\section{The Sigma Global Optimization Package}

	The Sigma algorithm \cite{Aluffi-Pentini:1988:ASE} finds the
global minimum of a function $f(\vec{x})$ of
$\vec{x}=(x_1,x_2,...,x_n)$ degrees of freedom.  The algorithm
generates a fixed number of trajectories over the $\vec{x}$
phase-space, with each trajectory updated by a move along the local
gradient of some random direction, plus an additional Gaussian noise
contribution for each of the $n$ components for $\vec{x}$.  In
addition, the package described here allows the user to specify a {\em
mask\/} which indicates which of the $n$ dimensions have known values
and should be treated as constants.

	From the user point of view, the basic unit of time in the
algorithm is the {\em observation period}, during which the set of
trajectories are computed.  At the end of each observation period one
of the trajectories is rejected, and another one selected for
bifurcation.  Also at the end of each observation period a set of
stopping criteria are examined to determine how the algorithm should
proceed.  The end of a series of observation periods constitute a {\em
trial}, and are deemed successful or unsuccessful, depending on
whether a solution is found.  The user may specify the algorithm to
execute until a specified number of successful trials are completed;
this is a possible insurance against a found minimum to be only a
local minimum.

	Three parameters determine the algorithm's macroscopic time
behaviour, and they are concerned with controlling the number of
observation periods.  These are:

\begin{itemize}
	\item The minimum number of observation periods to run for.
The author has observed that setting this quantity to be too low meant
the algorithm was unable to achieve `equilibrium', and local minima
were found too readily.

	\item The maximum number of observation periods to run for.
For large dimensional cases, or where the optimisation function is
hard, this value has to be set to be large.  On the other hand, it may
be more useful to have the algorithm `give up' and try again with
another starting condition.

	\item The number if observation periods to increment for, once
the minimum number has been reached.

\end{itemize}

	The following terms are used by Sigma:

\begin{itemize}

	\item {\em Uniform Stop}, refers to when all the trajectories
in a given trial return the same (within numerical accuracy) minimum
value.  Note that for degenerate minima, the {\em minimiser\/} for
each of the trajectories may not be all equal.

	\item {\em Successful Stop}, refers to when a uniform stop has
yielded the minimum value (numerically) equal to the {\em best\/}
minimum found since the algorithm started.

\end{itemize}


	The variance of the maximum likelihood function for two
parameters $\theta_a$ and $\theta_b$ is given by \cite{book:gpy:1984}
the expression:

\begin{equation}
V_{ab}^{-1} = - \left< \frac{\partial^2 \ln{\cal L}}
			{\partial\theta_a\partial\theta_b} \right>
\end{equation}

where $<{\cdots}>$ is the average over the liklihood function $\cal
L$, over all the observation variables.  This is equivalent to the
expression evaluated, at the {\em minimum\/} point\footnote{Since
that is what Sigma seeks to find.}:

\begin{equation}
\label{eq:Vab-1}
V_{ab}^{-1} = \left.\frac{\partial^2 \ln{\cal L}}
			{\partial\theta_a\partial\theta_b}
			\right|_{\hat{\theta}_a,\hat{\theta}_b}
\end{equation}

The matrix formed by the elements $V_{ab}^{-1}$ must then be inverted,
with the diagonal elements of the inverse matrix giving the variance
of a given parameter.  In the code, the elements~(\ref{eq:Vab-1}) are
computed by central finite difference, and the resulting
positive-definite matrix inverted by Cholesky decomposition
\cite{book:isd-ame-jkr:1986}.


\bibliography{biblio}\bibliographystyle{alpha}



\end{project}
\end{document}


