# Finance Models¶

## Brownian Motions: Quick Introduction¶

The law of motion for stocks is often based on a geometric Brownian motion, i.e.,

$$dS(t) = \mu S(t) \; dt + \sigma S(t) \; dB(t), \quad S(0)=S_0.$$

This is a stochastic differential equation (SDE), because it describes random movement of the stock $S(t)$. The coefficient $\mu$ determines the drift of the process, and $\sigma$ determines its volatility. Randomness is injected by Brownian motion $B(t)$. This is more general than a deterministic differential equation that is only a function of time, as with a bank account, whose accretion is based on the equation $dy(t) = r \;y(t) \;dt$, where $r$ is the risk-free rate of interest.

The solution to a SDE is not a deterministic function but a random function. In this case, the solution for time interval $h$ is known to be

$$S(t+h) = S(t) \exp \left[\left(\mu-\frac{1}{2}\sigma^2 \right) h + \sigma B(h) \right]$$

The presence of $B(h) \sim N(0,h)$ in the solution makes the function random. The presence of the exponential return makes the stock price lognormal. (Note: if random variable $x$ is normal, then $e^x$ is lognormal.)

Re-arranging, the stock return is

$$R(t+h) = \ln\left( \frac{S(t+h)}{S(t)}\right) \sim N\left[ \left(\mu-\frac{1}{2}\sigma^2 \right) h, \sigma^2 h \right]$$

Using properties of the lognormal distribution, the conditional mean of the stock price becomes

$$E[S(t+h) | S(t)] = e^{\mu h}$$

The following R code computes the annualized volatility $\sigma$.

The parameter $\mu$ is also easily estimated as

So the additional term $\frac{1}{2}\sigma^2$ does matter substantially.

## Monte Carlo Simulation¶

It is easy to simulate a path of stock prices using a discrete form of the solution to the Geometric Brownian motion SDE.

$$S(t+h) = S(t) \exp \left[\left(\mu-\frac{1}{2}\sigma^2 \right) h + \sigma \cdot e\; \sqrt{h} \right]$$

Note that we replaced $B(h)$ with $e \sqrt{h}$, where $e \sim N(0,1)$. Both $B(h)$ and $e \sqrt{h}$ have mean zero and variance $h$. Knowing $S(t)$, we can simulate $S(t+h)$ by drawing $e$ from a standard normal distribution.

## Vectorization¶

The same logic may be used to generate multiple paths of stock prices, in a vectorized way as follows. In the following example we generate 3 paths. Because of the vectorization, the run time does not increase linearly with the number of paths, and in fact, hardly increases at all.

The plot code shows how to change the style of the path and its color.

If you generate many more paths, how can you find the probability of the stock ending up below a defined price? Can you do this directly from the discrete version of the Geometric Brownian motion process above?

## Bivariate random variables¶

To convert two independent random variables $(e_1,e_2) \sim N(0,1)$ into two correlated random variables $(x_1,x_2)$ with correlation $\rho$, use the following trannsformations.

$$x_1 = e_1, \quad \quad x_2 = \rho \cdot x_1 + \sqrt{1-\rho^2} \cdot x_2$$

We see that these are uncorrelated random variables. Now we create a pair of correlated variates using the formula above.

Algebraic and numerical check

Check algebraically that $E[x_i]=0, i=1,2$, $Var[x_i]=1, i=1,2$. Also check that $Cov[x_1,x_2]=\rho = Corr[x_1,x_2]$.

We know that $E[e_i]=0, i=1,2$, $Var[e_i]=1, i=1,2$. Also

• $E(e_i^2) = Var(e_i) + E(e_i)^2 = 1 + 0 = 1$.
• $E(e_1 \cdot e_2) = Cov(e_1,e_2) + E(e_1)E(e_2) = 0 + 0 = 0$.
• $E(x_1) = E(e_1) = 0$.
• $E(x_2) = E(\rho \cdot e_1 + \sqrt{1-\rho^2}\cdot e_2)$.
• $Var(x_1) = E(x_1^2)-E(x_1)^2 = E(e_1^2) - 0 = 1-0 = 1$.
• $Var(x_2) = E(x_2^2)-E(x_2)^2 = E((\rho \cdot e_1 + \sqrt{1-\rho^2}\cdot e_2)^2)-0 = \rho^2 E(e_1^2)+(1-\rho^2)E(e_2^2) = 1$
• $Cov(x_1,x_2) = E(e_1 e_2) - E(e_1)E(e_2) = E[e_1(\rho \cdot e_1 + \sqrt{1-\rho^2}\cdot e_2)] = \rho E(e_1^2) + (1-\rho^2)E(e_1 \cdot e_2) = \rho$

## Multivariate random variables¶

These are generated using Cholesky decomposition. We may write a covariance matrix in decomposed form, i.e., $\Sigma = L\; L'$, where $L$ is a lower triangular matrix. Alternatively we might have an upper triangular decomposition, where $U= L'$. Think of each component of the decomposition as a square-root of the covariance matrix.

We now compute the Cholesky decomposition of the covariance matrix.

The Cholesky decomposition is now used to generate multivariate random variables with the correct correlation.

In the last calculation, we confirmed that the simulated data has the same covariance matrix as the one that we generated correlated random variables from.

## Portfolio Computations¶

Let's enter a sample mean vector and covariance matrix and then using some sample weights, we will perform some basic matrix computations for portfolios to illustrate the use of R.

The expected return of the portfolio is

And the standard deviation of the portfolio is

We thus generated the expected return and risk of the portfolio, i.e., the values $0.078$ and $0.187$, respectively.

We are interested in the risk of a portfolio, often measured by its variance. As we had seen in a previous chapter, as we increase $n$, the number of securities in the portfolio, the variance keeps dropping, and asymptotes to a level equal to the average covariance of all the assets. It is interesting to see what happens as $n$ increases through a very simple function in R that returns the standard deviation of the portfolio.

We now apply it to increasingly diversified portfolios.

We can plot this to see the classic systematic risk plot. Systematic risk declines as the number of stocks in the portfolio increases.

## Optimal Portfolio¶

We will review the notation one more time. Assume that the risk free asset has return $r_f$. And we have $n$ risky assets, with mean returns $\mu_i, i=1...n$. We need to invest in optimal weights $w_i$ in each asset. Let $w = [w_1, \ldots, w_n]'$ be a column vector of portfolio weights. We define $\mu = [\mu_1, \ldots, \mu_n]'$ be the column vector of mean returns on each asset, and ${\bf 1} = [1,\ldots,1]'$ be a column vector of ones. Hence, the expected return on the portfolio will be

$$E(R_p) = (1-w'{\bf 1}) r_f + w'\mu$$

The variance of return on the portfolio will be

$$Var(R_p) = w' \Sigma w$$

where $\Sigma$ is the covariance matrix of returns on the portfolio. The objective function is a trade-off between return and risk, with $\beta$ modulating the balance between risk and return:

$$U(R_p) = r_f + w'(\mu - r_f {\bf 1}) - \frac{\beta}{2} w'\Sigma w$$

The f.o.c. becomes a system of equations now (not a single equation), since we differentiate by an entire vector $w$:

$$\frac{dU}{dw'} = \mu - r_f {\bf 1} - \beta \Sigma w = {\bf 0}$$

where the RHS is a vector of zeros of dimension $n$. Solving we have

$$w = \frac{1}{\beta} \Sigma^{-1} (\mu - r_f {\bf 1})$$

Therefore, allocation to the risky assets

• Increases when the relative return to it $(\mu - r_f {\bf 1})$ increases.
• Decreases when risk aversion increases.
• Decreases when riskiness of the assets increases as proxied for by $\Sigma$.

Example

## What if we reduced beta?¶

Notice that the weights in stocks scales linearly with $\beta$. The relative proportions of the stocks themselves remains constant. Hence, $\beta$ modulates the proportions invested in a risk-free asset and a stock portfolio, in which stock proportions remain same. It is as if the stock versus bond decision can be taken separately from the decision about the composition of the stock portfolio. This is known as the two-fund separation property, i.e., first determine the proportions in the bond fund vs stock fund and the allocation within each fund can be handled subsequently.

## Interest-rate processes (O-U process)¶

Here are the equations for the Ornstein-Uhlenbeck (O-U) process, for mean-reverting interest rates, $r(t)$:

$$dr(t) = k(\theta -r(t)) dt + \sigma dZ(t), \quad r(0)=r_0$$

where $k$ is the rate of mean reversion, $\theta$ is the long run mean of the process, and $\sigma$ is its volatility parameter. The solution to this stochastic differential equation is:

$$r(t) = r(0) e^{-kt} + \theta (1-e^{-kt}) + \sigma \int_0^t e^{-k(t-s)}dZ(s)$$

Note that the solution is a random function. We can use the solution to determine the conditional moments of the process as follows:

$$E[r(t) | r(0)] = r(0) e^{-kt} + \theta (1-e^{-kt})$$$$Var[r(t)|r(0)] = \frac{\sigma^2}{2b}(1-e^{-2kt})$$

For simulation purposes, we discretize the solution as follows:

$$r(t) = r(0) e^{-kt} + \theta (1-e^{-kt}) + \sigma e^{-0.5kt}\cdot z\sqrt{t}$$

## Cox-Ingersoll-Ross (CIR) model¶

The O-U process permits negative interest rates. The CIR model handles this issue with a small modification:

$$dr(t) = k(\theta -r(t)) dt + \sigma \sqrt{r} \cdot dZ(t), \quad r(0)=r_0$$$$r(t) = r(0) e^{-kt} + \theta (1-e^{-kt}) + \sigma e^{-0.5kt}\cdot z\sqrt{t\cdot r(t)}$$

## Estimating the historical parameters for interest rates¶

For the O-U process, a simple regression may be used to estimate its parameters. Discretizing the process leads to the following equation:

$$r(t+h) - r(t) = k(\theta - r(t)) h + \sigma z \sqrt{h}, \quad z \sim N(0,1)$$

We write $y(t)=r(t+h) - r(t)$, $\alpha = k \theta h$, $\beta = -kh$, and $e(t) = \sigma z(t) \sqrt{h}$. Then the equation above is written as

$$y(t) = \alpha + \beta \cdot r(t) + e(t)$$

which is a regression equation from which we can recover the coefficients of the model as follows:

\begin{eqnarray} k &=& -\frac{\beta}{h} \\ \theta &=& \frac{\alpha}{kh} = -\frac{\alpha}{\beta} \\ Var[e(t)] &=& \sigma^2 h \end{eqnarray}

Therefore, also

$$\sigma = \sqrt{\frac{1}{h}\cdot Var[e(t)]} = \frac{1}{\sqrt{h}}\cdot \sigma_e$$

Run the estimation.

## Antithetic Variance Reduction¶

The antithetic variable for a normal random variable $x$ is $x' = -x$. For a uniform random variable $u$, the antithetic is $u' = (1-u)$. In both cases, we are assured that the covariance between a variable and its antithetic is $Cov(x,x') \leq 0$. We can repeat the simulation model using the antithetic approach. The simulated estimates will have lower variance.

We see that the variance has come down with the antithetic approach.

## Control Variate Technique (CVT)¶

This is an approach where a Monte Carlo pricing model $P$ is improved by benchmarking off a similar pricing model $Q$ for which a closed-form solution is known. The idea is to discover the simulation bias in model $Q$ and use it to correct the bias in model $P$.

Suppose the Monte Carlo estimate of the pricing problem is denoted $P$ and that of the analogous problem is denoted $Q$. The closed-form solution for the latter model is denoted $Q_0$. Error from simulation is proportional to $(Q-Q_0)$. We use this error to correct the price $P$ to a new price $P'$ as follows:

$$P' = P + (Q_0 - Q)$$

For this approach to work it must be that $Var(P') \leq Var(P)$, i.e., we achieve a reduction in pricing error. For this to be the case, the following condiiton must be satisfied:

$$Var(P') = Var(P + (Q_0 - Q)) = Var(P-Q) = Var(P) + Var(Q) - 2 Cov(P,Q)$$

This implies that for $Var(P')<Var(P)$ it must be that

$$Cov(P,Q) > \frac{1}{2}Var(Q)$$

In other words, model $Q$ must be closely related to model $P$ for this approach to be effective in reduction of pricing error.

A good example in option pricing for this approach is in the pricing of Asian options (i.e., options on the average stock price). There is no closed-form option pricing model for options on the arithmetic average price (model $P$) but there is for options on the geometric average price (model $Q$). The latter may be used to improve the pricing of the former.

## American Options: Simulation with Optimization¶

American options are usually priced using backward recursion, whereas simulation runs forward in time. How might we adapt forward simulation to price American put options?

We note that an American put is exercised early at time $t$ when the stock price $S$ drops low enough that the gains from exercise $K-S(t)$ are greater than the expected future gains from exercise. The locus of points $B(t)$ where current gains from exercise are equal to future gains from exercise is known as the "early-exercise boundary" of the American put option. If we are given the boundary, then it is easy to simulate the stock price and every time it hits the boundary, we exercise and take profits. The average of the present value of these profits over many simulated paths is the price of the American put option.

But, we do not have the early-exercise boundary! So, in order to use simulation, we search over many boundaries till we find the one that gives the highest option value. We posit the following functional form for the boundary: $B(T-t) = K \exp(-\alpha (T-t))$. This simple functional form has a single parameter $\alpha$. By varying $\alpha$ we get a family of boundaries to search over. The program code is as follows.

Let's run this for different $\alpha$ to find the optimal early-exercise boundary and the price of the American put.

Tha value here is greater than that of the European put option.

## American Options: Binomial Trees¶

We can reverify these results using a binomial tree approach to pricing the American put.

Run the program to get the correct American put price.

As we see, the price is close to the one obtained by simulation.

## Longstaff-Scwartz Least Squares Model¶

In an important paper, Francis Longstaff and Eduardo Schwartz developed a simple and interesting approach to pricing American options using regressions. We present the example from their paper here, see @doi:10.1093/rfs/14.1.113.