%Paper: ewp-fin/9702003
%From: jvic@tiac.net
%Date: Tue, 18 Feb 97 01:19:06 CST
%Date (revised): Thu, 31 Jul 97 01:25:15 CDT


\documentstyle[12pt]{article}

\title{Discount-Bond Derivatives on\\
a Recombining Binomial Tree$^1$}
%\thanks{Draft \#2. Sections 3 and 4
%Neumann boundary conditions has been added to this version.}}


\author{J. Chalupa \\ 
Box 82 \\
Princeton, MA 01541 USA \\
jvic@tiac.net}
 \date{}
\maketitle
\begin{abstract}
\noindent Interest-rate derivative models governed by 
parabolic partial differential equations (PDEs) are studied 
with discrete-time recombining binomial trees. For the 
B\"uhler-K\"asler discount-bond model, the expiration value 
of the bond is a limit point of tree sites. Representative 
calculations give a close approximation to the continuum results. 
Next, situations are considered in which spatial 
inhomogeneity of the drift velocity can cause binomial 
jump probabilities to become negative. 
When the continuous-time boundary 
conditions are applied near the tree points at which this
occurs, good agreement 
is obtained with with Hull and White's 
explicit-finite-difference treatment of the Cox-Ingersoll-Ross model. 
Finally, to mimic the effect of a drift-velocity divergence 
which prevents interest rates from becoming negative, 
Neumann boundary conditions are applied in the Vasicek model. 
Discrete-time computations are performed for a mean-reverting 
situation and for a case with constant negative 
short-rate drift; the ensuing bond values have nonnegative 
interest rates and forward rates. The results are compared 
with the Vasicek solution and with the leading term 
in a spectral expansion.  
 
\end{abstract}


\begin{document}
\footnotetext[1]{Draft No. 2. Sections 3 and 4 on 
Neumann boundary conditions are new. Additional references 
have been included, and some material from the previous 
draft has been moved to the Appendix.}
\newpage
{\noindent
\bf\large 1. Introduction\\}

This working paper applies binomial-tree models to calculate 
interest-rate 
derivative valuations. Particular attention is paid to the 
interplay between boundary conditions and the approach to the 
continuous-time limit. A natural variable for modelling a stock 
of price $S$ is $\ln S$, which is unbounded ($-\infty < \ln S 
< \infty$). Since interest rates are in principle 
nonnegative\footnote{\noindent In 
practice, models permitting negative interest rates should 
produce a negligible probability for such rates. For situations 
in which the predicted probability is not negligible, it 
is appropriate to bar negative rates outright.}, 
a discount-bond 
price B has an upper bound and $\ln B$ is confined to s 
half-space. The continuous-time valuations considered 
here are described by parabolic partial differential equations 
(PDEs); the correct behavior near the half-space 
boundary must be built into 
a tree model's passage to the continuous-time limit.  

The paper is organized as follows. Section 2 considers the 
B\"uhler-K\"asler discount-bond model\cite{buhlerkasler}. 
The binomial tree has 
a limit point which makes the boundary unattainable in a 
finite number of jumps. The Cox-Ingersoll-Ross (CIR) 
model\cite{cir} 
is treated in Section 3. When CIR is mapped onto a model with 
constant variance, Neumann boundary conditions (vanishing 
slope) hold in the new variable. Although these boundary conditions 
are not explicitly stressed in Nelson and Ramaswamy's binomial 
analysis\cite{nelson}, they can be applied naturally in the 
binomial-tree formulation; it is argued 
that the Hull-White modification of the explicit 
finite-difference method\cite{finitedifferencehull} converges to them. 
In Section 4, Neumann conditions are applied to the Vasicek 
model\cite{vasicek} as a way of 
precluding negative interest rates. The interest-rate-dependent 
mean-reversion term is suppressed, and a computed example 
has nonnegative discount-bond rates for constant {\it negative} 
drift velocity. Section 5 presents concluding 
remarks. The tree structure for a nonconstant variance 
has been derived concisely by use 
of Ito's lemma\cite{finitedifferencehull,nelson}, the Appendix 
gives alternative derivations, one of which admits a 
simple multifactor generalization\cite{jcmultifactor}.

\newpage  


{\noindent
\bf\large 2. Inaccessible Boundary: B\"uhler-K\"asler\\}

In the "direct" approach to bond valuation, the bond price 
is the fundamental random variable. 
Rady and Sandmann's article\cite{radysandmann} 
reviews this line of effort, especially key contributions 
by Ball and Torous\cite{balltorous}, by Briys, Crouhy and Sch\"obel
\cite{briysetal}, and by 
B\"uhler and K\"asler\cite{buhlerkasler}. The simplest 
lattice for a discrete-time formulation is a recombining 
binomial tree. 
The resulting picture has the same relationship 
to the continuous-time model as the Cox-Ross-Rubinstein 
binomial-tree description\cite{coxrossrubinstein} has to 
Black-Scholes\cite{blackscholes}.

Rady and Sandmann\cite{radysandmann} caution that 
in bond-price direct models, an 
absorbing boundary and nonzero probability 
of vanishing interest rate can arise. On a binomial tree, such 
behavior could be precluded if the boundary were a limit point 
which is inaccessible in a finite number of steps. Figure 1 
depicts such a tree, on which computations can readily be 
performed.   
Representative numerical results will be presented for the 
B\"uhler-K\"asler model, whose continuum limit  has an 
economic and probabilistic "clean bill of health". Moreover, 
the tree structure can be determined analytically. 



\begin{figure*}[bhtp]
\begin{center}
\unitlength=1mm
\begin{picture}(150,100)(-15,0)
\multiput(0,0)(0,10){10}{
\multiput(0,0)(10,0){5}{\circle*{1}}  }
\multiput(49.9,0)(0,10){10}{\circle*{1} }
\multiput(59.8,0)(0,10){10}{\circle*{1} }
\multiput(69.5,0)(0,10){10}{\circle*{1} }
\multiput(78.7,0)(0,10){10}{\circle*{1} }
\multiput(86.9,0)(0,10){10}{\circle*{1} }
\multiput(93.1,0)(0,10){10}{\circle*{1} }
\multiput(96.87,0)(0,10){10}{\circle*{1} }
\multiput(99.51,0)(0,10){10}{\circle*{1} }
\multiput(99.82,0)(0,10){10}{\circle*{1} }
\multiput(0,0)(0,10){9}{
\multiput(5,5)(10,0){5}{\circle*{1}}  }
\multiput(54.9,5)(0,10){9}{\circle*{1}}
\multiput(64.7,5)(0,10){9}{\circle*{1}}
\multiput(74.2,5)(0,10){9}{\circle*{1}}
\multiput(83.0,5)(0,10){9}{\circle*{1}}
\multiput(90.26,5)(0,10){9}{\circle*{1}}
\multiput(95.26,5)(0,10){9}{\circle*{1}}
\multiput(99.21,5)(0,10){9}{\circle*{1}}
\multiput(99.95,5)(0,10){9}{\circle*{1}}
\end{picture}
\end{center}
\caption{{\bf Schematic diagram of a binomial tree.} 
Time steps occur in the vertical direction, and the horizontal 
direction depicts price. 
Hops occur between nearest neighbors on adjacent time steps. The 
depicted tree becomes uniform to the left and has a limit point 
on the right. Thus there is an infinite number of nodes both 
to the left and to the right of any given node. }
\end{figure*}



B\"uhler and K\"asler value bond derivatives by solving the PDE

\begin{equation}
\frac{\partial f}{\partial t} + 
\frac{1}{2}\sigma^2\frac{\partial^2 f}{\partial z^2}-
\frac{1}{2}\sigma^2\frac{\partial f}{\partial z} = 0  
\label{blsch}
\end{equation} 

\noindent with variance 

\begin{equation}
\sigma^2(z)=\sigma_0^2 (1-e^z)^2.
\end{equation}

\noindent Equation (\ref{blsch}) is the Black-Scholes 
equation with a price-dependent variance. It is solved in the 
half-space $z \le 0$ with the boundary condition 

\begin{equation}
f(z=0,t)=f(z=0,T).
\end{equation}

\noindent Derivative securities assume their expiration value 
at $z=0$.


The tree structure is found by inverting

\begin{equation}
\xi=n\sqrt{\tau} = \int_{z_0}^z \frac{dz'}{\sigma(z')} 
\label{solution}
\end{equation}

\noindent for the z-coordinate of the $n$th 
node\cite{nelson,finitedifferencehull}. For the 
B\"uhler-K\"asler model, the solution is

\begin{equation}
z=-\ln\left[1+e^{-n\sigma_0\sqrt{\tau}}\left(e^{-z_0}-1
\right)\right].
\end{equation}

\noindent For $n\rightarrow-\infty$, the nodes become evenly 
spaced, but for $n\rightarrow+\infty$, $z=0$ is an accumulation 
point of nodes. In random-interest-rate models like the 
Dothan\cite{dothan} and Black-Karasinski\cite{blackkarasinski} 
models, 
an accumulation point would arise at vanishing short-term rate.) 
The accumulation point would not exist if $1/\sigma(z)$ were 
integrable near $z=0$; examples of such situations will 
be examined in the next two sections. For constant $\sigma$, 
Rady and Sandmann show that the boundary conditions on the PDE 
can be interpreted in terms of an absorbing boundary at $z=0$ 
and a finite probability for a vanishing interest rate; these are 
absent for the B\"uhler-K\"asler model. The absence of an 
absorbing barrier is consistent with the limit point at $z=0$. 

Table 1 shows representative numerical results showing good 
agreement with the B\"uhler-K\"asler call value
\cite{radysandmann,buhlerkasler}
\begin{eqnarray}
f(z,t) &=& (1-e^\epsilon) e^z 
N\left[\frac{1}{\sigma_0\sqrt{T-t}}\left(
\ln\frac{e^z (1-e^\epsilon)}{(1-e^z)e^\epsilon} 
  +\frac{\sigma_0^2 (T-t)}{2}\right)\right] \nonumber \\
&~~& -e^\epsilon (1-e^z) 
N\left[\frac{1}{\sigma_0\sqrt{T-t}}\left(
\ln\frac{e^z (1-e^\epsilon)}{(1-e^z)e^\epsilon} 
  -\frac{\sigma_0^2 (T-t)}{2}\right)\right]
\end{eqnarray}

\noindent where the strike price is $e^\epsilon$ and 
$T$ is the expiration time. For the given 
parameter values, the at-the-money B\"uhler-K\"asler call 
is more than a factor of $5$ less than the Black-Scholes call. 
The difference is due to the price dependence of the variance. 
This price dependence is captured by the nonuniform binomial 
tree. 

\begin{table} 
\begin{center}
\begin{tabular}{|l|c|c|}\hline
~~ & \shortstack{$$z=-.200000\\(n=0)$$} & 
\shortstack{$$z=-.214697\\(n=-10)$$} \\ \hline
$T-t=50\tau$ & .995 & .965 \\ \hline
$T-t=100\tau$ & .9975 & .990 \\ \hline
$T-t=200\tau$ & .999 & .997 \\ \hline
\end{tabular}
\end{center}
\caption{
{\bf Ratio of binomial calls to the B\"uhler-K\"asler value.}
The $\ln$ of the strike price is 
$\epsilon=-\frac{1}{5}$, and the node index $n$ 
is zero at $z_0=-\frac{1}{5}$. The standard deviation per time step 
is $\sigma_0 = \frac{3}{20 \cdot 365^{1/2}}$, and the time step is 
$\tau=1$.}
\end{table}
~~\\


{\noindent
\bf\large 3. Neumann Boundary Conditions: Cox-Ingersoll-Ross\\}

It is known that artificial negative jump probabilities can arise 
in literal discrete-time treatments of models with large 
drift velocities. The Cox-Ingersoll-Ross (CIR) model provides an 
example. CIR is a one-factor model driven by a random short-term 
interest $r$ whose variance is proportional to $r$. For vanishing 
price of risk and a drift velocity $v=a(b-r)$ with $a,b > 0$, the 
valuation equation for a derivative security is\cite{cir} 

\begin{equation}
\frac{\partial f}{\partial t} +
\frac{\sigma^2}{2}r\frac{\partial^2 f}{\partial r^2} + 
a(b-r)\frac{\partial f}{\partial r} - r f = 0
\label{cir1}
\end{equation}

\noindent The derivative $f$ can be priced in terms of 
(a presumably complete set of) eigenfunctions 
$\psi$ generated by the substitution 
$f \rightarrow e^{-\lambda (T-t)} \psi(r,\lambda)$. The 
transformation of variables $\xi=2\sqrt{r}/\sigma$ produces 

\begin{equation}
\frac{1}{2}\frac{\partial^2 \psi}{\partial \xi^2} +
\left[\left(\frac{2 a b}{\sigma^2}-
\frac{1}{2}\right)\frac{1}{\xi}-\frac{1}{2} a \xi\right] 
       \frac{\partial \psi}{\partial \xi} + 
\left(-\frac{1}{4}\sigma^2 \xi^2  +\lambda \right)\psi 
= 0 . 
\label{cir2}
\end{equation}

\noindent Values of $\lambda$ are determined by the boundary 
conditions. These are that $\psi$ is regular at $r=0$ and 
vanishes at $r=\infty$. Examination of the power series 
solution\cite{mf} 
indicates that the term of order $\xi$ is absent. 
Thus, the regularity requirement at $r=0$ becomes the 
Neumann boundary condition 

\begin{equation}
\frac{\partial \psi}{\partial \xi}=0
\end{equation}

\noindent Since $r\propto \xi^2$, this is satisfied by the 
CIR discount-bond value $A(t)e^{-B(t) r}$. 

How is the Neumann boundary condition to be implemented on the 
binomial tree associated with the time-dependent version of 
equation \ref{cir2}? The nodes are at 
$\xi=n\sqrt{\tau}\quad (n=0,1,2,...)$. For $n=0$ and perhaps some 
neighboring sites, and for sufficiently large values of $n$, one 
hopping probability is negative (and the other is 
greater than one). A natural procedure is to apply the boundary 
condition at the values $n^*$ of $n$ at which the probabilistic 
interpretation breaks down\footnote{By no means is this the only 
choice. For example, $v$ could be regularized, e.g. by 
$\frac{v\sqrt{\tau}}{\sigma}\rightarrow 
\frac{v\sqrt{\tau}}{\sqrt{\sigma^2+v^2\tau}}$.}; 
using the condition $\partial f(\xi=0,t)/\partial \xi=0$ rather 
than $f(\xi=0,t)=0$ permits $f$ to remain nonzero at finite $r$. 
Table 2 shows a comparison 
with the exact CIR results and Hull and White's 
explicit-finite-difference solution\cite{finitedifferencehull}. 
The Neumann conditions are given the form 
$f(n^*\sigma\sqrt{\tau},t)=f((n^*+1)\sigma\sqrt{\tau},t)$. 
This leads 
to two interleaving binomial trees which interact at their 
boundaries. The time step $\tau$ was chosen identical 
to Hull and White's, and the $\xi$-step size is 
$1/\sqrt{3}$ smaller than that on their trinomial tree; this is 
partially offset by the fact that 
computing a trinomial jump requires fewer multiplications than 
computing a binomial one. 

\begin{table} 
\begin{center}
\begin{tabular}{|l|c|c|c|c|c|}\hline
~~&r=.06&r=.08&r=.10&r=.12&r=.14 \\ \hline
~~&.6631&.6353&.6086&.5830&.5585 \\
$T-t=5$ yrs&(.6631)&(.6353)&(.6086)&(.5830)&(.5585) \\
~~&((.6624))&((.6345))&((.6078))&((.5823))&((.5578)) \\ \hline

~~&.4091&.3898&.3714&.3538&.3371 \\
$T-t=10$ yrs&(.4092)&(.3898)&(.3714)&(.3538)&(.3371) \\
~~&((.4083))&((.3889))&((.3705))&((.3529))&((.3362)) \\ \hline

~~&.2502&.2382&.2268&.2159&.2055 \\
$T-t=15$ yrs&(.2502)&(.2382)&(.2268)&(.2159)&(.2055) \\
~~&((.2494))&((.2374))&((.2259))&((.2151))&((.2147)) \\ \hline

~~&.1528&.1454&.1385&.1318&.1255 \\
$T-t=20$ yrs&(.1528)&(.1454)&(.1385)&(.1318)&(.1255) \\
~~&((.1521))&((.1448))&((.1378))&((.1311))&((.1248)) \\ \hline
\end{tabular}
\end{center}
\caption{
{\bf CIR Results.} Hull and White's parameters are 
adopted: $a=.4$, $b=.1$, 
$\sigma=.06$, and $\tau=.05$. The top entry in each box is the 
present calculation; linear interpolation was used on the 
tree values. The second (in parentheses) is Hull-White's 
explicit-finite-difference result, and the third 
((in double parentheses)) is the exact CIR result quoted 
by Hull and White.}
\end{table}


Actually, the Neumann boundary condition at $r=0$ is virtually 
implicit in the binomial formulation of CIR. 
Consider the continuous-time 
limit $\tau \rightarrow 0$; in this limit, 
$v(\xi =n\sqrt{\tau})\sqrt{\tau}$ approaches a constant 
for fixed $n$ 
and so do the corresponding jump probabilities. Ignoring the 
$1/(1+r\tau)$ term which is of higher order in $\xi$, one sees that 
for adjacent tree sites
\begin{eqnarray}
f(n\sqrt{\tau},t) &=& p(n) f((n+1)\sqrt{\tau},t+\tau)\nonumber\\
&~&\qquad +(1-p(n)), 
\end{eqnarray} 

\noindent where $p(n)$ is the jump probability at site $n$. 
The point is that $f(r\propto\xi^2,t)$ is assumed regular and 
smooth near $r=0$, but this cannot be true for $p(n)$ because 
of the $1/\xi$ term in the drift velocity. Expanding about 
$\xi=0$ to leading order in $\sqrt{\tau}$, one finds

\begin{equation}
\frac{\partial f(\xi=0,t)}{\partial \xi} \simeq 
(2 p(n)-1) \frac{\partial f(\xi=0,t+\tau)}{\partial \xi}.
\end{equation}

\noindent For positive drift velocity, $2 p -1 > 0$ holds, and, 
for $p(n)<1$, under repeated time steps a finite 
$\partial f(\xi \sim \sqrt{\tau},t)/\partial\xi$ is rapidly driven 
to zero. 

Similar considerations hold for the Hull-White 
explicit-finite-difference treatment of CIR. For example, 
consider the relation
\begin{eqnarray}
f(n\Delta\phi,t)&=&p_n^+ f((n+1)\Delta\phi,t+\tau)\nonumber\\
&~&\quad+(1-p_n^+ - p_n^-) f(n\Delta\phi,t+\tau)\nonumber\\
&~&\qquad +p_n^- f((n-1)\Delta\phi,t+\tau),
\end{eqnarray}

\noindent in which the $p_n^\pm$'s are the probabilities of 
up or down jumps at site $n$ and $\Delta\phi$ 
is the $\xi$-step size. Since $\Delta\phi=O(\sqrt{\tau})$, 
for $\tau\rightarrow 0$ at fixed $n$, this equation can be 
expanded about $\xi=0$. At first order in $\Delta\phi$, there 
ensues 
\begin{eqnarray}
n\Delta\phi \frac{\partial f(\xi=0,t)}{\partial \xi} &\simeq& 
(n+p_n^+ - p_n^-) \Delta\phi 
\frac{\partial f((n+1)\Delta\phi,t+\tau)}{\partial \xi}(0,t+\tau) \\
&\simeq& 
(n+p_n^+ - p_n^-) \Delta\phi 
\frac{\partial f(\xi=0,t)}{\partial \xi}.
\label{hwneumann}
\end{eqnarray}

\noindent Since $p_n^\pm$ approaches a constant in the 
limit considered, Neumann boundary conditions hold if 
the Taylor expansion is 
valid\footnote{Equation (\ref{hwneumann}) leaves open the 
possibility that $\partial f(\xi \rightarrow 0,t)/\partial\xi$ 
grows under iteration. For the functional form of Hull and 
White's jump probabilities and their bounds on the variables in 
the functions, it can be shown that this is not the case 
at $n=n^*$.}---not surprisingly, given 
the quantitative success of Hull-White reproduced 
in Table 2.\\

{\noindent
\bf\large 4. Neumann Boundary Conditions: Vasicek\\}

The power-series solution\cite{mf} of equations of the form 
(\ref{cir2}) indicates that the Neumann boundary condition 
describes the solution which is regular at the origin 
when the drift velocity has a pole there. In effect, a pole 
with positive residue prevents $\xi$ from becoming negative. 
It is natural to consider such a divergent drift 
velocity in a constant-variance interest rate model, i.e. in 
the Vasicek model, whose variants and 
generalizations\cite{holee,hullwhitestuff,blackkarasinski} 
are used to infer derivative valuations from empirical 
term structure data. The robustness of the results is of 
particular interest: when the residue of the pole is tuned 
toward zero, the limiting case is a "regular" Vasicek model 
with Neumann boundary conditions. Also, for constant $\sigma$, 
if Neumann conditions are adopted for the forward 
equation $\frac{\partial f}{\partial t} + 
\frac{1}{2}\sigma^2\frac{\partial^2 f}{\partial r^2} +
v(r) \frac{\partial f}{\partial r} = 0$, the adjoint boundary 
condition is $\frac{1}{2}\frac{\partial f}{\partial r} - v f =0$, 
which indicates the conservation of probability---time-independent 
normalization---in the backward (adjoint) equation 
$-\frac{\partial f}{\partial t} + 
\frac{1}{2}\sigma^2\frac{\partial^2 f}{\partial r^2} -
\frac{\partial v(r) f}{\partial r} = 0$. 

This section 
presents binomial-tree computations on the Vasicek 
model; the results will be 
interpreted with the continuous-time theory. First, a 
standard scenario with strong mean reversion will be 
considered; then, to highlight the role of boundary 
conditions, a situation without mean reversion and with 
{\it negative} interest-rate drift will be analyzed. 
The intention (at least in the present draft) is to 
demonstrate the procedure rather than 
to present a survey of the $(\sigma,b,a,\tau)$ 
parameter space. 

The valuation equation is 

\begin{equation}
\frac{\partial f}{\partial t} + 
\frac{1}{2}\sigma^2\frac{\partial^2 f}{\partial r^2} +
a(b-r) \frac{\partial f}{\partial r} -rf = 0.
\label{vasicekeqn}
\end{equation}

\noindent This equation is usually analyzed in the 
domain $-\infty < r < \infty$, where its behavior 
simplifies\cite{vasicek,jamshidian}. 
The mean reversion implied by the convention 
$a,b>0$ typically keeps the contribution of negative 
rates small, but this need not be the case for sufficiently 
volatile assets or for nonstandard values of $a$ and $b$. 
The solution can be expressed as an 
eigenfunction expansion:

\begin{equation}
f(r,t)=\sum_{n=0}^\infty \psi(r,\lambda_n) 
e^{-\lambda_n (T-t)}.
\end{equation}

\noindent The eigenfunction terms are the solutions of 
equation (\ref{vasicekeqn}) which vanish at $r=\infty$; 
the eigenvalues are determined by the boundary condition at 
$r=0$. The $\psi_\lambda$'s can be expressed in terms of 
confluent hypergeometric (parabolic cylinder) functions or Airy 
functions\cite{abramowitz}. 
Because (\ref{vasicekeqn}) is not a Sturm-Liouville 
equation, a weight factor must be used to orthonormalize the 
$\psi_\lambda$'s. Interest rates will be negative if the smallest 
$\lambda$ appearing in the expansion is negative. 


\begin{table} 
\begin{center}
\begin{tabular}{|l|c|c|c|c|c|c|}\hline
~~&r=.04&r=.06&r=.08&r=.10&r=.12&r=.14 \\ \hline
~~&.8048&.7555&.7093&.6658&.6251&.5868 \\
$T-t=5$ yrs&(.8163)&(.7390)&(.6688)&(.6052)&(.5477)&(.4956) \\
~~&((.8047))&((.7554))&((.7091))&((.6657))&((.6249))&((5866)) 
\\ \hline

~~&.6365&.5839&.5356&.4913&.4507&.4134 \\
$T-t=10$ yrs&(.6395)&(.5789)&(.5239)&(.4742)&(.4291)&(.3883) \\
~~&((.6364))&((.5836))&((.5353))&((.4910))&((.4503))&((.4130)) 
\\ \hline

~~&.5003&.4552&.4141&.3766&.3425&3115 \\
$T-t=15$ yrs&(.5010)&(.4535)&(.4105)&(.3715)&(.3361)&((.3042)) \\
~~&((.5003))&((.4550))&((.4134))&((.3762))&((.3421))&((.3111)) 
\\ \hline

~~&.3925&.3560&.3229&.2928&.2654&.2407 \\
$T-t=20$ yrs&(.3925)&(.3553)&(.3216)&(.2910)&(.2633)&(.2383) \\
~~&((.3926))&((.3559))&((.3226))&((.2924))&((.2651))&((.2403)) 
\\ \hline
\end{tabular}
\end{center}
\caption{
{\bf Results for Mean-Reverting Vasicek.} 
The following parameters are adopted: $a=.2$, $b=.05$, 
$\sigma=.01$, and $\tau=.04$. The top entry in each box is the 
binomial calculation; linear interpolation was used on the 
tree values. The second (in parentheses) is the one-mode 
approximation, and the third ((in double parentheses)) is the 
Vasicek solution for $-\infty < r < \infty$.}
\end{table}


Table 3 presents a canonical case with strong mean reversion.
Neumann conditions are imposed at $r=0$ and at the value of 
$r$ at which a jump probability becomes negative. The numerical 
results are in good agreement with the Vasicek solution; the 
dynamics is driven by the mean reversion and boundary effects 
are minor. The results 
are also consistent with a one-term truncation of the 
eigenfunction expansion; as expected, the agreement improves 
when time to maturity increases and the leading transient 
dominates. 


\begin{table} 
\begin{center}
\begin{tabular}{|l|c|c|c|c|c|c|}\hline
~~&r=.04&r=.06&r=.08&r=.10&r=.12&r=.14 \\ \hline
~~&.9059&.8356&.7600&.6883&.6231&.5640 \\
$T-t=5$ yrs&(.9039)&(.8202)&(.7168)&(.6036)&(.4902)&(.3840) \\
~~&((.9297))&((.8412))&((.7612))&((.6887))&((.6232))&((.5639)) 
\\ \hline

~~&.8736&.7974&.7007&.5978&.4992&.4120 \\
$T-t=10$ yrs&(.8817)&(.8001)&(.6992)&(.5888)&(.4781)&(3746) \\
~~&((1.1237))&((.9200))&((.7537))&((.6167))&((.5049))&((.4134)) 
\\ \hline

~~&.8505&.7732&.6771&.5718&.4666&.3691 \\
$T-t=15$ yrs&(.8600)&(.7804)&(.6820)&(.5743)&(.4637)&(.3654) \\ 
~~&((1.7883))&((1.3248))&((.9814))&((.7271))&((.5386))&((.3990)) 
\\ \hline

~~&.8254&.7504&.6570&.5544&.4511&.3543 \\
$T-t=20$ yrs&(.8389)&(.7612)&(.6652)&(.5602)&(.4549)&(.3564) \\
~~&((3.7934))&((2.5430))&((1.7046))&((1.1426))&((.7659))&((.5134)) 
\\ \hline
\end{tabular}
\end{center}
\caption{
{\bf Vasicek Results for Constant 
Drift.} The following parameters are 
adopted: $a=0$, $v=a b=-.01$, 
$\sigma=.01$, and $\tau=.04$. The top entry in each box is the 
binomial calculation; linear interpolation was used on the 
tree values. The second (in parentheses) is the one-mode 
approximation, and the third ((in double parentheses)) is the 
Vasicek solution for $-\infty < r < \infty$.}
\end{table}



Next, a nonstandard situation will be considered in which the 
boundary condition is the only factor working against negative 
interest rates. Equation (\ref{vasicekeqn}) will be considered 
in the limit $(a\rightarrow 0, a b\rightarrow v)$, and 
the interest-rate drift $v$ will be taken as negative. 
Negative interest rates are not apparent in 
Table 4. In the set of eigenfunctions  

\begin{equation}
\psi(r,\lambda) = e^{-r v/\sigma^2} Ai\left(
\frac{r-\lambda+\frac{v^2}{2 \sigma^2}}
{(\frac{1}{2}\sigma^2)^{1/3}}\right), 
\end{equation}

\noindent the asymptotic properties of the Airy function in the limit 
$\frac{v}{\sigma^{4/3}}\rightarrow -\infty$, lead to the estimate
$\lambda\approx \frac{\sigma^2}{2 |v|}$, in good agreement with the 
Newton's-method result $\lambda\simeq .00497$. (For $v=0$, the 
leading $\lambda$ 
is proportional to the Airy function's largest root, which is 
negative; for $v\rightarrow+\infty$, $\frac{v^2}{2\sigma^2}-
\lambda$ is proportional to the 
zero of $Ai'(x)$, which is also negative.)  \\

The absence of negative eigenvalues for constant drift raises the 
issue for $v=a(b-r)$ and more general cases.
For example, when $a>0$, the appropriate 
eigenfunction is 

\begin{eqnarray}
f(r,\lambda)&=&\pi e^{-r/a}
\left[\frac{M\left(c,\frac{1}{2},\frac{a}{\sigma^2}
  (r+\frac{\sigma^2}{a^2}-b)^2\right)}
    {\Gamma(\frac{1}{2}+c)\Gamma(\frac{1}{2})}\right. \nonumber\\
&~&\quad-\left.\frac{\sqrt{a}}{\sigma}(r+\frac{\sigma^2}{a^2}-b)
\frac{M\left(\frac{1}{2}+c,\frac{3}{2},\frac{a}{\sigma^2}
   (r+\frac{\sigma^2}{a^2}-b)^2\right)}
         {\Gamma(c)\Gamma(\frac{3}{2})}\right],
\end{eqnarray}

\noindent where

\begin{equation}
c=-\frac{1}{2 a}\left[\frac{\sigma^2}{2 a^2}-b+\lambda\right]. 
\end{equation}

\noindent and $f\rightarrow 0$ as $r\rightarrow \infty$; a similar 
solution exists for $a<0$. The leading 
eigenvalue $\lambda=b-\frac{\sigma^2}{2 a^2}$ in 
the Vasicek solution 
satisfies the equation $c=0$, which can be shown to 
determine an eigenvalue, i.e. $\frac{1}{\Gamma(c)}=0$, 
in the limit $\sigma\rightarrow 0$. For 
$b=\frac{\sigma^2}{a^2}$, examination of the $\Gamma$-functions 
indicates that the solutions 
to $\frac{\partial f(r=0,\lambda)}{\partial r} =0$ must have 
negative $c$ and positive $\lambda$. 
At present, the spectrum associated with an arbitrary value of 
$(a,b)$ must be determined on an individual basis, although so 
far all prospective negative roots of 
$\frac{\partial f(r=0,\lambda)}{\partial r} =0$ have turned out to 
be weakly positive or asymptotes. 

Nevertheless, a qualititive argument about the spectrum can be made. 
Consider the forward equation 
$-\frac{\partial f}{\partial t} + 
\frac{1}{2}\sigma^2\frac{\partial^2 f}{\partial r^2} -
\frac{\partial v(r) f}{\partial r} = 0$ 
mentioned at the beginning of this section. Since the (adjoint) 
boundary conditions give conservation of probability, the equation 
is assumed to have well-behaved solutions with a nonnegative 
spectrum. Therefore, the backward equation  

\begin{equation}
\frac{\partial f}{\partial t} + 
\frac{1}{2}\sigma^2\frac{\partial^2 f}{\partial r^2} +
v(r) \frac{\partial f}{\partial r} = 0,
\label{adjointprobability}
\end{equation}

\noindent with Neuman boundary conditions, is also well-behaved.

Consider a binomial-tree treatment of this equation. If 
the initial value of the solution is nonnegative, the properly 
bounded jump probabilities assure that it will be nonnegative at 
all prior times. Introducing the damping factor 
$\frac{1}{1+ r\tau} \leq 1$ at each jump in 
equation (\ref{adjointprobability}) produces a 
binomial-tree formulation of 

\begin{equation}
\frac{\partial f}{\partial t} + 
\frac{1}{2}\sigma^2\frac{\partial^2 f}{\partial r^2} +
v(r) \frac{\partial f}{\partial r} - r f = 0, 
\label{theequation}
\end{equation}

\noindent which is the equation of interest. If the temporal 
divergences associated with negative eigenvalues are not present 
in (\ref{adjointprobability}), they will not arise in 
(\ref{theequation}). Note that this rationale breaks down if 
$r$ is allowed to be negative. 

~\\
{\noindent
\bf\large 5. Summary and Conclusion\\}

Three things have been done in this paper. First, the 
discrete-time B\"uhler=K\"asler discont-bond model has 
been studied on a binomial tree which turned out to have 
a limit point. Second, it was noted that when the 
Cox-Ingersoll-Ross model is mapped onto a constant-variance 
model, the appropriate boundary condition in the new price 
variable is the Neumann condition; binomial-tree calculations 
are consistent with, and clarify the success of, the 
Hull-White explicit-finite-difference method. Third, as a 
heuristic way of modelling a short-range divergent drift 
which keeps interest rates nonnegative, Neumann conditions 
were imposed in the Vasicek model. In a typical mean-reverting  
situation, the effects are minor. When mean reversion is 
turned off and a constant negative interest-rate drift 
imposed, the consequences in a sample case are dramatic: 
the discount-bond valuation appears well-conditioned, 
whereas the Vasicek solution is dominated by 
negative rates. If mean reversion is not essential 
to obtain {\it prima facie} acceptable results from Vasicek, 
the model and its variants may accommodate a larger 
category of fitting functions, and be applicable to 
securities with greater volatility, than would be the 
case otherwise. 


\newpage
\noindent{\bf\Large Appendix\\}

The tree structure has been concisely 
derived\cite{nelson,finitedifferencehull} by use 
Ito's lemma. This Appendix derives it by formulating the 
appropriate continuous-time limit of a binomial process. 
A check using the PDE commences after equation (\ref{drat}); 
it is the starting point for a simple generalization to the 
multifactor case. 

In the binomial processes of interest, the price of 
interest-rate derivatives depends on the value of a 
random variable $z$ which can be, for 
instance, an interest rate or a bond price. Changes in $z$ occur 
in discrete time steps $\tau$.
At time $t$ and $z=z_0$, $z$ 
changes to $z_0+a_+$ with probability $p$ and to 
$z_0-a_-$ with probability $1-p$. 
Derivative securities can be 
valued with price-of-risk methods\cite{hull,merton}. It is 
convenient to adapt the terminology of Ref. \cite{jc}.
The value of a derivative 
security $f(z,t)$ is determined from its values one time step 
closer to expiration by 

\begin{equation}
\frac{
\left\langle f(z+\Delta z,t+\tau \right\rangle_{\Delta z} - 
(1+R)f(z,t) }{{\cal R}[f]} = \lambda\tau^{1/2}, 
\label{valuationequation}
\end{equation}

\noindent where $\lambda$ is the price of risk and $R$ is 
the interest. 
In principle, these quantities can depend on 
both $z$ and $t$. The expectation value is 

\begin{equation}
\left\langle f(z_0)\right\rangle \equiv p f(z+a_+) + 
(1-p)f(z-a_-)\qquad\qquad (a_+ + a_- > 0),  
\end{equation}

\noindent and the risk metric ${\cal R}$ is taken as 
quadratic:
\begin{eqnarray}
\left\vert{\cal R}[f]\right\vert
&=& \left\langle (f(z+\Delta z)-
\left\langle f \right\rangle)^2 \right\rangle^{1/2} \\
&=& \sqrt{p(1-p)}\left\vert f(z+a_+)-f(z-a_-)\right\vert
\end{eqnarray}

\noindent The Black-Scholes ${\cal R}$ has the sign of 
$\partial f/\partial z$, and in the present situation, 
${\cal R}$ is  assigned the sign of $f(z+a_+) - f(z-a_-)$. 
The solution to equation (\ref{valuationequation}) is 
\samepage{\begin{eqnarray}
f(z,t) &=& \frac{1}{1+R}\left[p(1-
\lambda\sqrt{\frac{\tau(1-p)}{p}})
     f(z+a_+,t+\tau)) + \right. \nonumber\\
&~&\left.\qquad\qquad\qquad
(1-p)(1+\lambda\sqrt{\frac{\tau p}{1-p}})
        f(z-a_-,t+\tau)\right] \\
&\equiv&  \frac{1}{1+R}\left[\tilde{p}f(z+a_+,t+\tau)) + 
(1-\tilde{p})f(z-a_-,t+\tau)\right].
\label{probeff}
\end{eqnarray} }

\noindent The functional form (\ref{probeff}) remains valid 
irrespective of the sign of ${\cal R}$.
If the underlying security is a bond of price $B$ and the 
random variable $z$ is $z=\ln B$, $\tilde{p}$ can be 
eliminated by requiring that $e^z$ be a solution of 
the equation, which transforms to

\begin{equation}
f(z,t)=\frac{(1+R-e^{-a_-})f(z+a_+,t+\tau)-
(1+R-e^{a_+})f(z-a_-,t+\tau)}{(1+R)(e^{a_+}-e^{-a_-})}.
\end{equation}

\noindent For a market consisting of the bond B, its derivatives, 
and time-independent cash, the above expression further 
simplifies to

\begin{equation}
f(z,t)=\frac{(1-e^{-a_-})f(z+a_+,t+\tau)+
(e^{a_+}-1)f(z-a_-,t+\tau)}{(e^{a_+}-e^{-a_-})}.
\label{bukassolution}
\end{equation}

\noindent after the interest rate is dropped. Equation 
(\ref{bukassolution}) describes a binomial hopping process 
with 

\begin{equation}
\frac{1-e^{-a_-}}{e^{a_+}-e^{-a_-}},\qquad
\frac{e^{a_+}-1}{e^{a_+}-e^{-a_-}}
\end{equation}

\noindent as transition probabilities.

As noted in Section 2, B\"uhler and K\"asler value 
bond derivatives by solving the PDE

\begin{equation}
\frac{\partial f}{\partial t} + 
\frac{1}{2}\sigma^2\frac{\partial^2 f}{\partial z^2}-
\frac{1}{2}\sigma^2\frac{\partial f}{\partial z} = 0  
\end{equation} 

\noindent subject to 
\begin{eqnarray}
\sigma^2(z)&=&\sigma_0^2 (1-e^z)^2 \\
f(z=0,t)&=&f(z=0,T).
\end{eqnarray}

\noindent Derivative securities assume their expiration value 
at $z=0$.

Can a binomial process be constructed which reduces to 
B\"uhler-K\"asler in the $\tau \rightarrow 0$ limit? A recombining 
tree is sought. If $\sigma$ depends on price but not on time, 
so will the lattice geometry. Let the price-dependent 
nodes of the tree be indexed by an integer $n$. At a given time 
step, all values of $n$ are even or odd; hops take place from site 
$n$ to sites $n \pm 1$.
Accordingly, the quantities $a_\pm$ depend on the index n. The 
constraint 

\begin{equation}
a_+(n) = a_-(n+1).
\end{equation}

\noindent holds because the tree structure is recombining and 
time-independent. 

A convenient feature of the B\"uhler-K\"asler model is that the 
variance approaches a constant values for large negative $z$. 
The tree becomes uniform in this domain with $a_+(n)\simeq
a_-(n)\simeq a \simeq \sigma_0 \tau^{1/2}$ as $n \rightarrow 
-\infty$ and $\tau \rightarrow 0$. For small $\tau$, the 
replacement 

\begin{equation} 
a_\pm(n) = \alpha_\pm \tau^{1/2}
\label{timepricedecouple}
\end{equation}

\noindent is reasonable; the uniform-tree results suggest that the
hopping probabilities should be expanded to $O(\tau^{1/2)}$. 
Thus the drift velocity is 

\begin{equation}
v_{drift}\tau \simeq -\frac{1}{2}\sigma^2(z(n))\tau 
\simeq -\frac{1}{2}a_+(n)a_-(n) 
=-\frac{1}{2} a_+(n) a_+(n-1). 
\end{equation}

\noindent The first equality above comes from the PDE; the second 
is a small-$a_\pm$ expansion of the tree model. 
If the variance rather than the drift velocity is 
computed, to leading order in $\tau$ the same relationship,

\begin{equation}
\sigma^2(z(n))\tau = a_+(n) a_+(n-1),
\label{recursionequation}
\end{equation}

\noindent ensues. Specifying one of the $a$'s determines all the 
rest. 

Position on the tree can be described in terms of $(z,t)$ or 
$(n,t)$. A natural question is whether a continuum model can 
be developed in the $\tau\rightarrow 0 $ limit. An obvious 
variable to try is 

\begin{equation}
\xi \equiv n \sqrt{\tau}. 
\end{equation}

\noindent Summations over n are thus Riemann sums:

\begin{equation}
\sqrt{\tau}\sum_{n}\psi(n\sqrt{\tau}) \longrightarrow 
\int d\xi \psi(\xi).
\end{equation}

\noindent Equation ({\ref{recursionequation}) for the $a_+$'s can 
be simplified by noting that 

\begin{equation}
a_+(n+1) \simeq a_+(n) \equiv a(n)
\end{equation}

\noindent to leading order in $\tau$ 
because $a_+(\xi+\tau) \simeq a_+(\xi)$.
The $\alpha_\pm$'s defined 
previously can be written as 

\begin{equation}
\alpha_\pm \simeq \alpha(\xi). 
\end{equation}

\noindent If the index $n$ is set to zero at $z_0$, i.e. 
$n(z_0)=0$, a $z$-position on the tree is determined by the number 
of $a$-steps from $z_0$:

\begin{equation}
z(\xi)-z_0 \simeq \sum_{n=0}^{n(z)} a_+(n) \simeq
\int_0^\xi d\xi' \alpha(\xi'). 
\end{equation}

\noindent Accordingly, to leading order in $\tau$, equation 
(\ref{recursionequation}) can be written as 

\begin{equation}
\frac{dz}{d\xi}=\sigma(z), 
\end{equation}

\noindent and inverting

\begin{equation}
\xi=n\sqrt{\tau} = \int_{z_0}^z \frac{dz'}{\sigma(z')} 
\label{drat}
\end{equation}

\noindent yields the $z$-coordinate of the $n$th node.


The foregoing discussion has been intuitive, and a consistency 
check is desirable. In fact, there is another interpretation 
which naturally encompasses models other than B\"uhler-K\"asler. 
The starting point is the continuum form 
$\xi=\int \frac{dz}{\sigma}$ 
for $\xi$ given in 
equation (\ref{solution}). If its coefficients are 
time-independent, a PDE of the form 

\begin{equation}
\frac{\partial f}{\partial t} +
\frac{\sigma^2}{2}\frac{\partial^2 f}{\partial z^2} + 
v_{eff} \frac{\partial f}{\partial z} - r f = 0
\end{equation}

\noindent can be rewritten as 

\begin{equation}
\frac{\partial f}{\partial t} +
\frac{1}{2}\frac{\partial^2 f}{\partial \xi^2} + 
(\frac{v_{eff}}{\sigma}-
  \frac{1}{2\sigma}\frac{\partial\sigma}{\partial\xi}) 
\frac{\partial f}{\partial \xi} - r f = 0
\label{xiequation}
\end{equation}

\noindent When 
$\int \frac{dz}{\sigma}$ diverges at zero and infinity, the 
relevant binomial process associated with equation 
(\ref{xiequation}) involves price-dependent hopping 
of the form (\ref{probeff}) on a uniform $\xi$-tree. 
The internode spacing 
is $\sqrt{\tau}$, and the hopping probability corresponding to 
the drift velocity of $\xi$ is 
\begin{eqnarray}
2 \tilde{p}-1&=&(\frac{v_{eff}}{\sigma}-
\frac{1}{2\sigma}\frac{\partial\sigma}{\partial\xi})\tau^{1/2}\\
&=&(\frac{v_{eff}}{\sigma}-
\frac{1}{2}\frac{\partial\sigma}{\partial z})\tau^{1/2}.
\end{eqnarray}

\noindent Consistency of the  $z$-tree and 
$\xi$-tree must be demonstrated. 
To order $\sqrt{\tau}$, the spacings are related by 
\begin{eqnarray}
\tau^{1/2}&=&\int_z^{z+a_+(z)} \frac{dz'}{\sigma(z')} 
\simeq \frac{1}{\sigma(z)}\int_z^{z+a_+(z)} dz' 
\left(1 - \frac{(z'-z)}{\sigma(z)}\frac{d\sigma(z)}{dz}\right)\\
a_+(z) &\simeq& \sigma(z)\tau^{1/2}+
\frac{1}{2}\sigma(z)\tau \frac{d\sigma(z)}{dz}.
\end{eqnarray}

\noindent The corresponding value of $a_-$ is 

\begin{equation}
a_-(z) \simeq \sigma(z)\tau^{1/2}-
\frac{1}{2}\sigma(z)\tau \frac{d\sigma(z)}{dz}.
\end{equation}

\noindent Because the relationships
\begin{eqnarray}
\sigma(z)^2 \tau &\simeq& a_+(z) a_-(z) \\
v_{eff}(z)\tau &\simeq& \tilde{p}(z) a_+(z) - 
 (1-\tilde{p}(z)) a_-(z)
\end{eqnarray}

\noindent hold to $O(\tau)$, the two trees are consistent. 
Discrete B\"uhler-K\"asler 
is recovered for the appropriate choice of $\tilde{p}$.



\newpage
\begin{thebibliography}{99}

\bibitem{buhlerkasler} W. B\"uhler, Zeitsch. 
Betriebswirtschaftliche Forsch.  $\underline{\hbox{10}}$, 
851 (1988); W. B\"uhler and J. K\"asler, Universit\"at 
Dortmund Working Paper (1989). I have used Rady and 
Sandmann's description of these papers. 

\bibitem{cir}
J.C. Cox, J.E. Ingersoll, and S.A. Ross, Econometrica 
$\underline{\hbox{53}}$, 385 (1985). 

\bibitem{nelson}
D.B. Nelson and K. Ramaswamy, Review of Financial Studies 
$\underline{\hbox{3}}$, 393 (1990).

\bibitem{finitedifferencehull}
J. Hull and A. White, J. Fin. Quant. Analysis 
$\underline{\hbox{25}}$, 87 (1990).

\bibitem{vasicek} O. Vasicek, J. Financial Economics
$\underline{\hbox{5}}$, 177 (1977).

\bibitem{jcmultifactor}
J. Chalupa,  Economics Working Paper Archive
ewp-fin/9706001 (1996) {\tt <}http://econwpa.wustl.edu{\tt >}

\bibitem{radysandmann} S. Rady and K. Sandmann, University of 
Bonn Discussion Paper 
B-212 (March 1995), published in Rev. of Futures Markets 
$\underline{\hbox{13}}$, 461 (1994). 

\bibitem{balltorous} C.A. Ball and W.N. Torous, J. Fin. Quant. 
Analysis  $\underline{\hbox{18}}$, 517 (1983); S. Cheng, J. 
Economic Theory $\underline{\hbox{53}}$, 185 (1991); A.G.Z 
Kemna, J.F.J. de Munnik, and A.C.F. Vorst, Erasmus Universiteit 
Rotterdam Discussion Paper, 1989, as summarized by Rady and 
Sandmann. 

\bibitem{briysetal} E. Briys, M. Crouhy, and R. Sch\"obel, 
J. Finance $\underline{\hbox{46}}$, 1879, (1991). 

\bibitem{coxrossrubinstein}
J.C. Cox, S. Ross, and M. Rubinstein, J. Financial Econ. 
$\underline{\hbox{7}}$, 229 (1979).

\bibitem{blackscholes}
F. Black and M. Scholes, J. Political Economy, 
$\underline{\hbox{81}}$, 637 (1973).

\bibitem{dothan}  L.U. Dothan, 
J. Financial Econ., $\underline{\hbox{6}}$, 59 (1978). 

\bibitem{blackkarasinski} 
F. Black and P. Karasinski, Financial Analysts Journal 
July-August 1991, p. 52. Cf. F. Black, E. Derman, and W. Toy, 
Financial Analysts Journal January-February 1990, p. 33.

\bibitem{mf}
P.M. Morse and H. Feshbach, {\it Methods of 
Theoretical Physics}, McGraw-Hill, New York, 1953. 

\bibitem{holee}
T.S.Y. Ho and S.-B. Lee, J. Finance 
$\underline{\hbox{41}}$, 1011 (1986).

\newpage

\bibitem{hullwhitestuff}
J. Hull and A. White, Review of Financial Studies 
$\underline{\hbox{3}}$, 573 (1990); J. Fin. Quant. Analysis 
$\underline{\hbox{28}}$, 235 (1993); J. Derivatives 
$\underline{\hbox{2}}$, 7 (1994); 
{\it ibid.} $\underline{\hbox{2}}$, 37 (1994); 
{\it ibid.} $\underline{\hbox{3}}$, 26 (1996).

\bibitem{jamshidian}
F. Jamshidian, J. Finance $\underline{\hbox{44}}$, 
205 (1989).

\bibitem{abramowitz}
M. Abramowitz and I.A. Stegun, {\it Handbook of Mathematical 
Functions}, Dover, New York, 1968.

\bibitem{hull}
J.C. Hull, {\it Options, Futures, and Other 
Derivative Securities}, Prentice Hall, Englewood Cliffs, 1993.

\bibitem{merton}
R.C. Merton, {\it Continuous-Time Finance},
Blackwell, Cambridge, 1994. 

\bibitem{jc}
J. Chalupa,  Economics 
Working Paper Archive
ewp-fin/9607009 (1996)
{\tt <}http://econwpa.wustl.edu/{\tt >}
(the email address in this paper is 
outdated); and 
references therein, particularly 
J.P. Bouchaud, G. Iori, and D. Sornette, RISK 
$\underline{\hbox{9}}$, 61 (1996) and 
E. Aurell and K. \.Zyczkowski, Economics 
Working Paper Archive ewp-fin/9601001 (1996), 
submitted to J. Political Economy. 

\end{thebibliography}

\end{document}