"> Kyo's Blog --> Kyo’s Blog | MBA Brand Manager Data Scientist

Kyo's Blog

MBA | Brand Manager | Data Scientist

Posted in Business, Analyics and tagged Business Analyics

This is a modified version of the report for Business Analyics in Bath Full Time MBA Class of 2020.

SALES FORECASTING 2

Introduction

In this section, it is assumed a pizza parlor business plan in the Bath University has lunched. The vital issue of the plan is also done so to find one or two pizzas which might be appealed by the Bath University students based on the attributions such as type of topping and base. Hence, it will be conducted an experiment to attempt diverse pizzas for the target consumers including the students. However, there are some constrictions for the plan, one of which is the maximum budget of £100 to purchase the pizzas for this experiment. One of the major tools to find the customer preferences can be conjoint analysis which assumed their preferences as a combination of the attributions having a certain degree (Janssens et al., 2008). Therefore, it should be decided that which attributions might be the most relevant to the purpose that is to find the most preferred among the students. The consideration relevant to the attributions was started from the research suggesting 5 attributions which are the crust, topping, cheese type, cheese quality and price (Lilien and Rangaswamy, 2002). The crust has 3 levels which are thin, thick, pan. Topping’s levels are pineapple, veggie, chicken, turkey. The levels of cheese type are Romano, mixed, mozzarella. These of cheese quality are 50g, 75g, 100g. £4.99, £6.99, £8.99 are the level of price. As a result, the total possible combinations of attribute levels are 324 (= 3 x 4 x 3 x 3 x 3). All of them cannot be conducted due to the time and budget restrictions, however, some researchers argued the effective number of part worth parameters to be approximated from conjoint data, the number of which is less than the sum of all levels of each attribution (Wierenga, 2008). Thus, if at most 16 (= 3 + 4 + 3 + 3 + 3) candidates are selected precisely, these alternatives can be represented all of the possible combinations of attribute levels (ibid). The selection usually is done by SPSS from all of attribute levels and therefore, the data format to save 16 alternatives representing all possibilities should be suitable for the application. Furthermore, the data structure of that should be 16 rows with at least 7 columns which are 5 attributions, the primary key and the preference score. The most well-liked pizza should the maximum preference score (one hundred) and the minimum one (zero) should be rated to the least popular pizza. Hence. Each pizza must have the preference score from 0 to 100. Conjoint analysis assumes the preference score can be represented the weighted sum of the utilities for each attribution (Wierenga, 2008). This means after collecting the score from the customers being tested, the software (SPSS) can discover weights and utilities of each attributions. In addition, the key is required to identify each rows, however, it is automatically inserted (see the leftist column of figure 1) when the data is inputted in the data format and then also accordingly the column named ‘STATUS’ with ‘Design’ as the value is inserted as a result of using the function to generate the orthogonal design (Janssens et al., 2008).

Figure 1

Furthermore, although it should be tested 16 type of pizzas to collect the target customer preference as the score, there is the budged limitation which the pizzas only can be purchased less than 100 pounds. This means if assumed the price of normal size pizza is 15 pounds, 6 pizzas are able to be bought for the experiment. Supposed the shop selected to buy the pizzas can allow to consume at least 3 type of the alternatives among 16 of that in the same pizza, 18 different type of pizza is available for the test. Therefore, all of alternatives are able to be examined with only 2 duplicates (16 = 18 - 2). Because of the fact that 6 students are applicable for the evaluation, if each student can eat and rate all of alternatives to cut each piece into 6 parts, the preferable scores of them are able to be collected. Figure 2 is assumed the result of the 16 alternative pizzas rating by 6 students (identified by id in the figure). For instance, the column named ‘sco1’ means the students (identified by id in the figure 2) rating score for the top row in the figure 1 (Card ID 1 row).

Figure 2

The utilities of each attributions are estimated by the software, SPSS. In Figure 3, it is seen apparently that the attribution ‘thin’ (assigned to ‘Crust’ level) has the highest utility (5.708) and the lowest utility (-4.010) is earned by the attribution ‘pineapple’ (belonged to ‘Topping’ level). Utilities

Figure 3

In addition, Figure 4 shows the averaged importance for each levels of the attributions which mean which levels can earn the highest utility. In the figure, the level ‘Topping’ has the highest one (22.503) and the lowest one (17.930) is gained by the level ‘Price’.

Figure 4

Furthermore, Figure 5 indicates that both the Pearson’s R (.864) and the Kendall’s tau (.639) which mean the correlations between the recognized and forecasted preferences are high (ibid). both correlations also are likely to be significant (respectively .000 and .000) because less than .05 significance means the null hypothesis of a null correlation is quite unlikely to be acceptable (ibid). Correlations2

Figure 5

Finally, Figures 6 through 10, the relative utility levels are demonstrated for each attribute and In Figure 11, the mean importance is displayed for each attribute (ibid). From Figures 6 through 11, it could be seen that the pizza parlor that offers their customers a thin of crust, Veggie of topping, Romano of cheese type, 50g of cheese quantity and pricing £8.99 might achieve the highest utility.

Figure 6 Figure 7 Figure 8 Figure 9 Figure 10 Figure 11

Reference

Posted in Business, Finance and tagged Business Finance

Organize and implement Hull and White (1994a) interest rate model

Derivation of Hull and White (1994a) interest rate model solutions and boundary conditions

  • Definition [Hull and White(1994a) interest rate model] \begin{eqnarray}\label{HullandWhite1994a} dr(t)=(\theta(t)-ar(t))dt+\sigma d\tilde{W}(t) ． \end{eqnarray}

However, $\tilde{W}(t)$ is Brownian motion under the risk-neutral measure $\mathbb{P}$. Also, $\theta (t)$, $a$, and $\sigma$ are positive-definite functions for time (Hull and White, 1994a).

  • Definition [Discount bond price evaluation formula] \begin{eqnarray}\label{DiscountBondPrice} B(t,T)=\tilde{\mathbb{E}}[e^{-\int_{t}^{T}R(s)ds}|{\cal F}(t)] ． \end{eqnarray}

Here, a discount bond is a contract in which $1$ is paid for a certain fixed maturity of $T$, and no payment is made before the maturity.

  • Definition [Yield from time $t$ to $T$] \begin{eqnarray} Y(t,T)=-\frac{1}{T-t}\log B(t,T)． \end{eqnarray}

This is equivalent to the following equation.

\begin{eqnarray} B(t,T)=e^{-Y(t,T)(T-t)}． \end{eqnarray}

Therefore, $B(T, T) =1$.

The random variable $e^{-\int_{t}^{T}R(s)ds}$ estimated in the equation \eqref{DiscountBondPrice} depends not only on $R(T)$ but also on the path. However, the path of $R$ to time $t$ depends only on $R(t)$ at time $t$. Therefore, since $R$ is a Markov process, the bond price $B(t, T)$ is considered to be a function of time $t$ and $R(t)$.

\begin{eqnarray} B(t,T)=f(t,R(t))． \end{eqnarray}

On the other hand, since the equation \eqref{DiscountBondPrice} is the expansion of $D(t)B(t,T)=\tilde{\Bbb{E}}[D(T)B(T,T))\lvert\cal{F}(t)]$, then $D(t)B(t,T)=D(t)f(t,T)$ is a martingale.

Therefore, in order for the no arbitrage condition to hold, the drift $(dt)$ term obtained by differentiating $D(t)f(t,T)$ must appear as $0$.

To find the derivative of $D(t)f(t,T)$, it is required first to find the discount process and its relational expression.

  • Definition [discount process] \begin{eqnarray} D(t)=e^{-\int_{0}^{t}R(s)ds}． \end{eqnarray}

Here, if $S(x)=e(-t)$ and $I(t)=\int_{0}^{t}R(s)ds$, then $D(t) = S(I(t))$. At this time, since $dI(t)=R(t)dt$, $dI(t)dI (t)=0$ and $S’(x)=-S (x)$, $S’‘(x) =S(x)$ then from Ito’s lemma,

\begin{eqnarray} dD(t)=dS(I(t)) =S’(I(t))dI(t)+\frac{1}{2}S’‘(I(t))dI(t)dI(t) =-S(I(t))dI(t) =-D(t)R(t)dt． \end{eqnarray}

therefore,

\begin{eqnarray} dD(t)=-D(t)R(t)dt． \end{eqnarray}

If the interest rate process $R(t)$ in finding the derivative of $D(t)f(t,T)$ follows the equation \eqref{HullandWhite1994a}, then From $dtdt=0$, $dtd\tilde{W}=0$, $d\tilde{W}d\tilde{W}=dt$ and $drdr= \sigma^2dt$

\begin{eqnarray} d(D(t)f(t,T)) =f(t,T)dD(t)+D(t)df(t,T) =f(t,T)(-D(t)R(t))+D(t){f_t(t,T)dt+f_r(t,T)dr+\frac{1}{2}f_rr(t,T)drdr} =D(t)(-rf(t,T)+f_t(t,T)+\frac{1}{2}\sigma^2 f_rr(t,T)+(\theta(t)-ar(t))f_r(t,T))dt+D(t)\sigma f_r(t,T)d\tilde{W}(t)． \end{eqnarray}

Therefore, $f(t, T))$ must be met.

\begin{eqnarray}\label{HullWhitePDE} f_t(t,T)+\frac{1}{2}\sigma^{2} f_rr(t,T)+(\theta(t)-ar(t))f_r(t,T)=rf(t,T) \end{eqnarray}

Also, since $1=B(T, T)=f(T, r)$, the terminal condition is obtained.

For all $r$, \begin{eqnarray} f(T,r)=1 \end{eqnarray}

Next, the partial differential equation of the equation \eqref{HullWhitePDE} is expanded to obtain an analytical solution. If you set $f_t(t, T)=e^{-rC (t, T) -A (t, T)}$,

\begin{eqnarray} f_t(t,r)=(-rC’(t,T)-A’(t,T))f(t,r)． f_r(t,r)=-C(t,T)f(t,r)． f_rr(t,r)=-C^{2}(t,T)f(t,r)． \end{eqnarray}

Where $C’(t,T)=\frac{\partial C(t,T)}{\partial t}$，$A’(t,T)=\frac{\partial A(t,T)}{\partial t}$． Therefore, the left side of the formula \eqref{HullWhitePDE} is

\begin{eqnarray} f_t(t,T)+\frac{1}{2}\sigma f_rr(t,T)+(\theta(t)-ar(t))f_r(t,T) =(-rC’(t,T)-A’(t,T)+\frac{1}{2}\sigma^{2}C^{2}(t,T)-(\theta(t)-ar(t))C(t,T))f(t,r)． \end{eqnarray}

Therefore, by transposing the right side to the left side and enclosing it with $r$, the following relational expression is obtained.

\begin{eqnarray} ((-C’(t,T)+aC(t,T)-1)r-A’(t,T)-\theta(t)C(t,T)+\frac{1}{2}\sigma^{2}C^{2}(t,T))f(t,r)=0． \end{eqnarray}

From the assumption, $f_t(t, T)=e^{-rC(t, T)-A(t, T)}>0$, and this relation must hold for all $r$, then this term for $r$ must be $0$.

\begin{eqnarray} C’(t,T)=aC(t,T)-1． A’(t,T)=-\theta(t)C(t,T)+\frac{1}{2}\sigma^{2}C^{2}(t,T)． \end{eqnarray}

Furthermore, $C(T, T)=A’(T, T)=A(T, T)=0$ is required to hold for the terminal condition.

A function that satisfies $C’(t, T)=aC(t, T)-1$, $C(T, T)=0$

\begin{eqnarray} C(t,T)=\int_{t}^{T}e^{-\int_{t}^{s}a(v)dv}ds． \end{eqnarray}

It is shown that this is a unique solution to be obtained.

\begin{eqnarray} \frac{d }{d s}(e^{-\int_{0}^{s}a(v)dv} C(s,T)) =e^{-\int_{0}^{s}a(v)dv}(-a(s)C(s,T)+C’(s,T) ) =e^{-\int_{0}^{s}a(v)dv}(-a(s)C(s,T)+aC(t,T)-1 ) =e^{-\int_{0}^{s}a(v)dv}． \end{eqnarray}

Integrating both sides of this equation from $t$ to $T$

\begin{eqnarray} e^{-\int_{0}^{T}a(v)dv}C(T,T)-e^{-\int_{0}^{t}a(v)dv}C(t,T)=-\int_{t}^{T}e^{-\int_{0}^{s}a(v)dv}ds \Leftrightarrow C(t,T)=e^{\int_{0}^{t}a(v)dv}\int_{t}^{T}e^{-\int_{0}^{s}a(v)dv}ds
=\int_{t}^{T}e^{-\int_{t}^{s}a(v)dv}ds． \end{eqnarray}

The expansion of this unique solution is as follows. However, $a(v)=a$ is set from the equation \eqref{HullandWhite1994a}.

\begin{eqnarray} C(t,T)=\int_{t}^{T}e^{-\int_{t}^{s}a(v)dv}ds =\int_{t}^{T}e^{-as+at}ds =\frac{1}{a}(1-e^{-a(T-t)})． \end{eqnarray}

Similarly, the function, that satisfies $A(T, T)=1$, $A’(t, T)=-\theta(t)C(t, T)+\frac{1}{2}\sigma^{2}C^{2}(t, T)$, is shown below as a unique solution.

\begin{eqnarray} A(t,T)=\int_{t}^{T}(\theta(t)C(t,T)-\frac{1}{2}\sigma^{2}C^{2}(t,T) )ds． \end{eqnarray}

Integrating both sides of this equation from $t$ to $T$ \begin{eqnarray} A(T,T)-A(t,T)=-\int_{t}^{T}\theta(t)C(t,T)ds+\frac{1}{2}\int_{t}^{T}\sigma^{2}C^{2}(t,T)ds \Leftrightarrow A(t,T)=\int_{t}^{T}(\theta(t)C(t,T)-\frac{1}{2}\sigma^{2}C^{2}(t,T))ds． \end{eqnarray}

The expansion of this unique solution is as follows. However, $a(v)=a$ is set from the equation \eqref{HullandWhite1994a}.

\begin{eqnarray} A(t,T)=\int_{t}^{T}(\theta(t)C(t,T)-\frac{1}{2}\sigma^{2}C^{2}(t,T) )ds =\int_{t}^{T}(\frac{\theta(t)}{a}(1-e^{-a(T-t)})-\frac{1}{2}\sigma^{2}\frac{1}{a^2}(1-e^{-a(T-t)} )^{2} )ds =\frac{\sigma^{2}}{2a^3}(-a(T-t)+2(1-e^{-a(T-t)})-\frac{1}{2}(1-e^{-2a(T-t)})) +\frac{1}{a}\int_{t}^{T}\theta(s)(1-e^{-a(T-s)})ds． \end{eqnarray}

Next, the parameter $\theta(t)$ in the equation \eqref{HullandWhite1994a} means the speed of mean reversion. It is required to match this parameter to the term structure of the quoted interest rate.

  • Definition [Instantaneous Forward Rate] Let $f^{M}(0, T)$ be the instantaneous forward rate of the market at maturity $T$ at time $0$. \begin{eqnarray} f^{M}(0,T)=-\frac{\partial \ln B^{M}(0,T)}{\partial T}． \end{eqnarray} Here, $B^{M}(0, T)$ is the discount bond price quoted in the market at maturity $T$.

Because it is $B(t,T)=f(t,r)=e^{-rC(t,T)-A(t,T)}$，

\begin{eqnarray} \ln B(0,T)=-rC(0,T)-A(0,T) =-\frac{r}{a}(1-e^{-aT} ) -\frac{\sigma^{2}}{2a^3}(-aT+2(1-e^{-aT})-\frac{1}{2}(1-e^{-2aT})) +\frac{1}{a}\int_{0}^{T}\theta(s)(1-e^{-a(T-s)})ds \Leftrightarrow \frac{\partial \ln B^{M}(0,T)}{\partial T}=re^{-aT}+\frac{\sigma^{2}}{2a^3}(-a+2ae^{-aT}+e^{-2aT})+\int_{t}^{T}\theta(s)e^{-a(T-s)}ds \Leftrightarrow f^{M}(0,T)=re^{-aT}+\frac{\sigma^{2}}{2a^3}(-a+2ae^{-aT}+e^{-2aT})+\int_{t}^{T}\theta(s)e^{-a(T-s)}ds． \end{eqnarray}

The following relationship is used to derive $\theta(T)$ from $f^{M}(0, T)$.

\begin{eqnarray} H(t)=\int_{0}^{t}h(s,t)dsarrow H’(t)=\int_{0}^{t}\frac{\partial h(s,t)}{\partial t}ds+h(t,t)． \end{eqnarray}

Here, $x(t)=re^{-at}+\int_{0}^{t}\theta(s)e^{-a(t-s)}ds$

\begin{eqnarray} x’(t)=-a(re^{-at}+\int_{0}^{t}\theta(s)e^{-a(t-s)}ds)+\theta(t) x’(t)=-ax(t)+\theta(t) \Leftrightarrow \theta(t)=x’(t)+ax(t) arrow \theta(T)=x’(T)+ax(T)． \end{eqnarray}

On the other hand, if it is set $g(t)=\frac{\sigma^2}{2a^2}(-1+2e^{-at}+e^{-2at})$

\begin{eqnarray} f^{M}(0,T)=x(T)+g(T) \Leftrightarrow x(T)=f^{M}(0,T)-g(T)． \end{eqnarray}

Therefore \begin{eqnarray} \theta(T)=\frac{\partial f^{M}(0,T)}{\partial T}-\frac{\partial g(T)}{\partial T}+ax(t)． \end{eqnarray}

Here, since it is $\frac{\partial g(T)}{\partial T}=-\frac{\sigma^2}{a}(e^{-aT}+e^{-2aT})$, the following equation is obtained after all.

\begin{eqnarray} \theta(T)=\frac{\partial f^{M}(0,T)}{\partial T}+\frac{\sigma^2}{a}(e^{-aT}+e^{-2aT})+af^{M}(0,T)-\frac{\sigma^2}{2a^2}(-1+2e^{-at}+e^{-2at})． \end{eqnarray}

Hull and White Tree formulation

In the following, to represent the interest rate model as the $3$ term tree, it is derived the relational expression between the parameters of the continuous time process and the parameters of the $3$ term process.

In the $3$ term process, it is expressed by $p_u$, $p_m$, and $p_d$, respectively, during the minute time $\Delta t$ as the interest rate rises by the spatial division width $\Delta x$, or stays at the same level, decreasing by $\Delta x$.

Furthermore, $\Delta t$ and $\Delta x$ cannot be selected independently and must satisfy the conditions for stability and convergence of the finite difference method $\Delta x=\sigma \sqrt{3\Delta t}$ (Hull and White, 1994).

Here, if the drift term of the formula \eqref{HullandWhite1994a} is set $\nu=\theta(t)-ar(t)$,

\begin{eqnarray}\label{averrageA} \mathbb{E}[\Delta x]=p_u(\Delta x)+p_m(0)+p_d(-\Delta x)=\nu\Delta t． \end{eqnarray} \begin{eqnarray}\label{varianceA} \mathbb{E}[\Delta x^2]=p_u(\Delta x^2)+p_m(0)+p_d(-\Delta x^2)=\sigma^2\Delta t+\nu^2\Delta t^2． \end{eqnarray} \begin{eqnarray}\label{SumProA} p_u+p_m+p_d=1． \end{eqnarray}

The expression \eqref{averrageA} and the expression \eqref{varianceA} are expanded as follows.

\begin{eqnarray}\label{averrageAA} p_u\Delta x-p_d\Delta x=\nu \Delta t． \end{eqnarray} \begin{eqnarray}\label{varianceAA} p_u\Delta x^2-p_d\Delta x^2=\sigma^2\Delta t+\nu^2\Delta t^2． \end{eqnarray}

Multiply the expression \eqref{averrageAA} by $\Delta x$ and add the expression \eqref{varianceAA},

\begin{eqnarray}\label{puA} 2p_u\Delta x^2=\nu^2\Delta t^2+\sigma^2\Delta t+\nu \Delta t\Delta x p_u=\frac{1}{2}(\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}+\frac{\nu \Delta t}{\Delta x} )． \end{eqnarray}

Also, substituting the expression \eqref{puA} into the expression \eqref{averrageAA},

\begin{eqnarray} p_u\Delta x-p_d\Delta x=\frac{1}{2}(\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}+\frac{\nu \Delta t}{\Delta x} )\Delta x-p_d\Delta x \Leftrightarrow p_d\Delta x=\frac{1}{2}(\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}-\frac{\nu \Delta t}{\Delta x} )\Delta x
p_d=\frac{1}{2}(\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}-\frac{\nu \Delta t}{\Delta x} )． \end{eqnarray}

And if it is substituted the expression \eqref{puA} and the expression \eqref{puA} into the expression \eqref{SumProA}

\begin{eqnarray}\label{pmA} p_u+p_m+p_d=\frac{1}{2}(\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}+\frac{\nu \Delta t}{\Delta x} )+p_m+\frac{1}{2}(\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}-\frac{\nu \Delta t}{\Delta x} ) =\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}+p_m \Leftrightarrow p_m=1-\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}． \end{eqnarray}

Hull And White (1994a) proposed a method for constructing an efficient and accurate short rate for the equation \eqref{HullandWhite1994a}. Also, it is assumed $\theta(t)=0$，$r(0)=0$.

\begin{eqnarray} dr(t)=ar(t)dt+\sigma d\tilde{W}(t) ． \end{eqnarray}

Therefore, note that $\nu=-ar$. In this simplified tree, the value of $r$ at the time $i\Delta t$ is $j\Delta t$. Since the tree is designed to reduce the amount of calculation, $j$ needs to be between $BottomNode[i]$ and $TopNode[i]$. The choice of which branching process to take at each node depends on the sign of $a$ and the size of $j$, and the range of change $\Delta x(=\Delta r)$ in the simplified $r$ process is chosen to match the expected value of $\mathbb{E}[\Delta x]$ and the variance $Var(\Delta x)$. Furthermore, when $a$ is positive, if $M$ is the expected value of the fluctuation range, when $j$ is the smallest integer larger than $-0.184/M$, the standard branching process is switched to the downward branching process. In order to maintain the symmetry, HW proposes to switch from the standard branching process to the upward branching process in the node $j$ with a negative sign for this node (Hull and White, 1994).

Here， \begin{eqnarray} \mathbb{E}[\Delta x]=\nu \Delta t =-ar\Delta t =Mr \Leftrightarrow M=-a\Delta t． \end{eqnarray}

Also， \begin{eqnarray} Var(\Delta x)=\mathbb{E}[\Delta x^2]-(\mathbb{E}[\Delta x])^2 =(\nu^2\Delta t^2+\sigma^2\Delta t)-(\nu \Delta t)^2 =\sigma^2\Delta t． \end{eqnarray} It becomes. Here,$M^2=a^2\Delta t^2$，$r=j\Delta r$,

\begin{eqnarray}\label{ExForPb1} \frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}=\frac{(-ar)^2\Delta t^2+\sigma^2\Delta t}{ (\sigma \sqrt(3\Delta t))^2} =\frac{(-a\Delta t)^2r^2+\sigma^2\Delta t}{ 3\sigma^2 \Delta t} =\frac{M^2(j\Delta r)^2+\sigma^2\Delta t}{ 3\sigma^2 \Delta t} =\frac{1}{3}\frac{j^2 M^2 \Delta r^2+\sigma^2\Delta t}{ \sigma^2 \Delta t} =\frac{1}{3}\frac{j^2 M^2 \Delta r^2}{ \sigma^2 \Delta t}+\frac{1}{3} =\frac{1}{3}\frac{j^2 M^2 (3\sigma^2 \Delta t)}{ \sigma^2 \Delta t}+\frac{1}{3} =j^2M^2+\frac{1}{3}． \end{eqnarray}

Therefore from the formula \eqref{ExForPb1}, \begin{eqnarray} p_m=1-\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2} =1-(j^2M^2+\frac{1}{3}) =\frac{2}{3}-j^2M^2． \end{eqnarray}

On the other hand, \begin{eqnarray}\label{ExForPb2} \frac{\nu\Delta t}{\Delta x}=\frac{(-ar)\Delta t}{\sigma\sqrt(3\Delta t)} =\frac{-Mr}{\sigma\sqrt(3\Delta t)} =\frac{-jr\Delta r}{\sigma\sqrt(3\Delta t)} =\frac{-jr(\sigma\sqrt(3\Delta t))}{\sigma\sqrt(3\Delta t)} =-jM． \end{eqnarray}

Therefore from the formula \eqref{ExForPb2},

\begin{eqnarray} p_u=\frac{1}{2}(\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}+\frac{\nu \Delta t}{\Delta x} ) =\frac{1}{2}( j^2M^2+\frac{1}{3} -jM ) =\frac{1}{6}+\frac{j^2M^2-jM}{2}． \end{eqnarray}

\begin{eqnarray} p_d=\frac{1}{2}(\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}-\frac{\nu \Delta t}{\Delta x}) =\frac{1}{2}(j^2M^2+\frac{1}{3}+jM) =\frac{1}{6}+\frac{j^2M^2+jM}{2}． \end{eqnarray}

Next, it is derived from the probabilities corresponding to the upward branching process.

\begin{eqnarray}\label{averrageB} \mathbb{E}[\Delta x]=p_u(2\Delta x)+p_m(\Delta x)+p_d(0)=\nu\Delta t． \end{eqnarray} \begin{eqnarray}\label{varianceB} \mathbb{E}[\Delta x^2]=p_u(2\Delta x)^2+p_m(\Delta x^2)+p_d(0)=\sigma^2\Delta t+\nu^2\Delta t^2． \end{eqnarray} \begin{eqnarray}\label{SumProB} p_u+p_m+p_d=1． \end{eqnarray}

The expression \eqref{averrageB} and the expression \eqref{varianceB} are expanded as follows.

\begin{eqnarray}\label{averrageBB} 2p_u\Delta x+p_m\Delta x=\nu \Delta t． \end{eqnarray} \begin{eqnarray}\label{varianceBB} 4p_u\Delta x^2+p_m\Delta x^2=\sigma^2\Delta t+\nu^2\Delta t^2． \end{eqnarray}

Subtracting the expression \eqref{averrageBB} from the expression \eqref{varianceBB} multiplied by $\Delta x$

\begin{eqnarray}\label{puB} 2p_u\Delta x^2=\nu^2\Delta t^2+\sigma^2\Delta t-\nu \Delta t\Delta x p_u=\frac{1}{2}(\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}-\frac{\nu \Delta t}{\Delta x} )． \end{eqnarray}

Also, if it is assigned the expression \eqref{puB} to the expression \eqref{averrageBB}, \begin{eqnarray}\label{pmB} 2p_u\Delta x+p_m\Delta x=\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}-\frac{\nu \Delta t}{\Delta x}+p_m\Delta x \Leftrightarrow p_m\Delta x=2\nu \Delta t-\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}
\Leftrightarrow p_m=\frac{2\nu \Delta t}{\Delta x}-\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}． \end{eqnarray}

And if it is substituted the expression \eqref{puB} and the expression \eqref{pmB} into the expression \eqref{SumProB} \begin{eqnarray}\label{pdB} p_u+p_m+p_d=\frac{1}{2}(\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}-\frac{\nu \Delta t}{\Delta x} )+\frac{2\nu \Delta t}{\Delta x}-\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}+p_d =\frac{3}{2}\frac{\nu \Delta t}{\Delta x}-\frac{1}{2}\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}+p_d \Leftrightarrow p_d=1-\frac{1}{2}( \frac{3\nu \Delta t}{\Delta x}-\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2} ) ． \end{eqnarray}

Therefore from the formula \eqref{ExForPb1}, \begin{eqnarray} p_u=\frac{1}{2}(\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}-\frac{\nu \Delta t}{\Delta x} ) =\frac{1}{2}( j^2M^2+\frac{1}{3}-(-jM)) =\frac{1}{6}-\frac{j^2M^2+jM}{2}． \end{eqnarray}

\begin{eqnarray} p_m=\frac{2\nu \Delta t}{\Delta x}-\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2} =-2jM-(j^2M^2+\frac{1}{3}) =\frac{1}{3}-j^2M^2-2jM． \end{eqnarray}

\begin{eqnarray} p_d=1-\frac{1}{2}( \frac{3\nu \Delta t}{\Delta x}-\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2} ) =1-frac{1}{2}( -3jM-j^2M^2-\frac{1}{3}) =\frac{7}{6}+\frac{j^2M^2+3jM}{2}． \end{eqnarray}

Then, the probability corresponding to the downward branching process is derived.

\begin{eqnarray}\label{averrageC} \mathbb{E}[\Delta x]=p_u(0)+p_m(-\Delta x)+p_d(-2\Delta x)=\nu\Delta t． \end{eqnarray} \begin{eqnarray}\label{varianceC} \mathbb{E}[\Delta x^2]=p_u(0)+p_m(-\Delta x)^2+p_d(-2\Delta x)^2=\sigma^2\Delta t+\nu^2\Delta t^2． \end{eqnarray} \begin{eqnarray}\label{SumProC} p_u+p_m+p_d=1． \end{eqnarray}

The expression \eqref{averrageC} and the expression \eqref{varianceC} are expanded as follows. \begin{eqnarray}\label{averrageCC} -p_m\Delta x-2p_d\Delta x=\nu \Delta t． \end{eqnarray} \begin{eqnarray}\label{varianceCC} p_m\Delta x^2+4p_d\Delta x^2=\sigma^2\Delta t+\nu^2\Delta t^2． \end{eqnarray}

Add the expression \eqref{varianceCC} multiplied by $\Delta x$ and the expression \eqref{averrageCC}, \begin{eqnarray}\label{pdC} 2p_d\Delta x^2=\nu^2\Delta t^2+\sigma^2\Delta t+\nu \Delta t\Delta x p_d=\frac{1}{2}(\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}+\frac{\nu \Delta t}{\Delta x} )． \end{eqnarray}

Also, if it is assigned the expression \eqref{pdC} to the expression \eqref{averrageCC}, \begin{eqnarray}\label{pmC} -p_m\Delta x-2p_d\Delta x=-p_m\Delta x-\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}-\frac{\nu \Delta t}{\Delta x} \Leftrightarrow -p_m\Delta x=2\nu \Delta t+\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}
\Leftrightarrow p_m=-\frac{2\nu \Delta t}{\Delta x}-\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}． \end{eqnarray}

And if it is substituted the expression \eqref{pdC} and the expression \eqref{pmC} into the expression \eqref{SumProC}, \begin{eqnarray}\label{puC} p_u+p_m+p_d=p_u-\frac{2\nu \Delta t}{\Delta x}-\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}+\frac{1}{2}(\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}+\frac{\nu \Delta t}{\Delta x} ) =p_u-\frac{3}{2}\frac{\nu \Delta t}{\Delta x}-\frac{1}{2}\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2} \Leftrightarrow p_u=1+\frac{3}{2}\frac{\nu \Delta t}{\Delta x}+\frac{1}{2}\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}． \end{eqnarray}

Therefore from the formula \eqref{ExForPb1}, \begin{eqnarray} p_d=\frac{1}{2}(\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2}+\frac{\nu \Delta t}{\Delta x} ) =\frac{1}{2}( j^2M^2+\frac{1}{3}-jM) =\frac{1}{6}-\frac{j^2M^2-jM}{2}． \end{eqnarray}

\begin{eqnarray} p_m=-\frac{2\nu \Delta t}{\Delta x}-\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2} =2jM-(j^2M^2+\frac{1}{3}) =-\frac{1}{3}-j^2M^2+2jM． \end{eqnarray}

\begin{eqnarray} p_u=1+\frac{3}{2}\frac{\nu \Delta t}{\Delta x}+\frac{1}{2}\frac{\nu^2\Delta t^2+\sigma^2\Delta t}{\Delta x^2} =1+\frac{3}{2}(-jM)+\frac{1}{2}(j^2M^2+\frac{1}{3}) =\frac{7}{6}+\frac{j^2M^2-3jM}{2}． \end{eqnarray}

Finally, the short rate process is introduced into the drift term $a_i$, which fluctuates with time. This drift term is chosen to be consistent with the discount bonds market quoted. \begin{eqnarray} P(0,i+1)=\int_{j=BottomNode[i]}^{TopNode[i]}Q_{i,j}e^{(-a_i+j\Delta r)\Delta t}． \end{eqnarray}

Here, $Q_{i,j}$ is a security that pays $1$ when the node $(i,j)$ is reached, and pays nothing at other times. Therefore, the drift term $ a_i $ is given as follows.

\begin{eqnarray} \ln P(0,i+1)=\ln \int_{j=BottomNode[i]}^{TopNode[i]}Q_{i,j}e^{j\Delta r\Delta t}-a_i\Delta t \Leftrightarrow a_i\Delta t=\ln \int_{j=BottomNode[i]}^{TopNode[i]}Q_{i,j}e^{j\Delta r\Delta t}-\ln P(0,i+1) \Leftrightarrow a_i=\frac{\ln \int_{j=BottomNode[i]}^{TopNode[i]}Q_{i,j}e^{j\Delta r\Delta t}-\ln P(0,i+1)}{\Delta t}． \end{eqnarray}

Hull and White Tree Implementation

It is given the simplified tree implementation mentioned above. The source code is written in matlab 2016b.

 
#container { float: left; margin: 0 -240px 0 0; width: 100%; } 
    for i=1:N-1 
        deltaX(i+1)=sigma*sqrt(3.*(t(i+1)-t(i)));
    end
    for i=1:N-1 
        for j=m(i):-1:-m(i) 
            M=-j.*deltaX(i).*a*(t(i+1)-t(i));
            V=sigma.^2.*(t(i+1)-t(i ));
            k(N+1-j, i)=round((j.*deltaX(i)+M)./deltaX(i+1)); 
            alpha=(j.*deltaX(i)+M-k(N+1-j, i).*deltaX(i+1))./deltaX(i+1);
            pu(N+1-j, i)=V./(2.*deltaX(i+1).^2)+(alpha.^2+alpha)./2;
            pd(N+1-j, i)=V./(2.*deltaX(i+1).^2)+(alpha.^2-alpha)./2; 
            pm(N+1-j, i)=1-V./(deltaX(i+1).^2)-alpha.^2;
        if j==max(m(i)) 
            m(i+1)=k(N+1-j,i)+1; 
        end
        end
    end
    for i=1:N-1
        for j=1:m(i)
            x(N+1+j, i)=x(N+1+j-1, i)-deltaX(i); 
            x(N+1-j, i)=x(N+1-j+1, i)+deltaX(i); 
        end
    end  
    for i=1:N-1 
        if i==2
            Q(N+2, i)=Q(N+1, i-1).*pu(N+1, i-1).*exp(-r(N+1, i-1).*(t(i)-t(i-1))); 
            Q(N+1, i)=Q(N+1, i-1).*pm(N+1, i-1).*exp(-r(N+1, i-1).*(t(i)-t(i-1))); 
            Q(N, i)=Q(N+1, i-1).*pd(N+1, i-1).*exp(-r(N+1, i-1).*(t(i)-t(i-1))); 
        elseif i>=3
            for j=m(i):-1:-m( i ) 
                for jj=m(i-1):-1:-m(i-1) 
                    Q(N+1+j, i)=Q(N+1+j, i)...
                    +Q(N+1+jj, i-1).*pu(N+1+jj, i-1).*exp(-r(N+1+jj, i-1).*(t(i)-t(i-1))).*(k(N+1+jj,i-1)+1==j)...
                    +Q(N+1+jj, i-1).*pm(N+1+jj, i-1).*exp(-r(N+1+jj, i-1).*(t(i)-t(i-1))).*(k(N+1+jj,i-1)==j)...
                    +Q(N+1+jj, i-1).*pd(N+1+jj, i-1).*exp(-r(N+1+jj, i-1).*(t(i)-t(i-1))).*(k(N+1+jj,i-1)-1==j);
                end 
            end 
        end
    if i==1
        g(i)=-log(P(i+1))./(t(i+1)-t(i)); 
    elseif i>=2 
        sum=0; 
        for j=m(i) :-1 :-m(i) 
            sum=sum+Q(N+1+j, i).*exp(-x(N+1+j, i).*(t(i+1)-t(i))); 
        end 
            g(i)=(log(sum)-log(P(i+1)))./(t(i+1)-t(i)); 
    end
        for j=m(i):-1:-m(i)
            r(N+1-j, i)=x(N+1-j, i)+g(i); 
            d(N+1-j, i)=exp(-r(N+1-j, i).*(t(i+1)-t(i))); 
        end 
    end

Reference

  • Hull, J. and A. White, 1994. “Numerical procedures for implementing term structure models I:Single-Factor Models,’’ Journal of Derivatives, 2, 1 (Fall 1994a) 7−16
  • Strickland, C. and Clewlow, L., 1998. Implementing Derivatives Models
  • Shreve, S., 2004. Stochastic Calculus for Finance II: Continuous-Time Models

Posted in Business, Finance and tagged Business Finance

Derivation of Black-Scholes PDE and its analytical solution by arbitrage pricing theory

In general, when the underlying asset price is given by a stochastic differential equation, the price of the derivative product is given by a partial differential equation. In this paper, it will be derived from the Black-Scholes PDE corresponding to the derivative (here, option) price valuation of the underlying asset (here, stock) modeled as the geometric Brownian motion. The valuation of an option is given by determining the initial funding needed to fully hedge the selling position of the option. Since the determined initial fund is the amount required to build a duplicate portfolio, it is called the no arbitrage price at time zero of the option. The price valuation method that duplicates options in stock market and money market account transactions is called arbitrage pricing theory (APT).

Suppose it is given a portfolio valued at $X(t)$ at any time $t$. Here, the composition of the given portfolio is as follows.

  • money-market account with interest rate $r$
  • A stock modeled by the geometric Brownian motion expressed by the following equation

\begin{eqnarray}\label{stockPrice} dS(t)=\alpha S(t)dt+\sigma S(t)dW(t)． \end{eqnarray}

Here, it is assumed that the investor who holds the given portfolio holds the stock in the unit of $\Delta(t)$ at any time $t$. At this time, it is assumed that the number of shares held $\Delta(t)$ is an ‘adapted stochastic process’. The remainder of the portfolio excluding the investment in stocks, $X(t)-\Delta(t)S(t)$, is invested in the money market account. Here, the small change $dX(t)$ of the portfolio value $X(t)$ at any time $t$ results in the following two factors.

  • Profit from holding shares $\Delta(t)dS(t)$
  • Profit from holding cash $r(X(t)-\Delta(t)S(t))dt$

That is, the small change in portfolio value $dX(t)$ is given by the following equation.

\begin{eqnarray}\label{portDiff} dX(t)=\Delta(t)dS(t)+r(X(t)-\Delta(t)S(t))dt． \end{eqnarray}

The equation \eqref{portDiff} is expanded as follows from the stochastic differential equation \eqref{stockPrice} of a given stock price. \begin{eqnarray}\label{portDiff2} dX(t)=\Delta(t)dS(t)+r(X(t)-\Delta(t)S(t))dt　 =\Delta(t)( \alpha S(t)dt+\sigma S(t)dW(t) )+r(X(t)-\Delta(t)S(t))dt =rX(t)dt+\Delta(t)( \alpha -r )S(t)dt+\Delta(t)\sigma S(t)dW(t)． \end{eqnarray}

Here, the three terms appearing in the small change $dX(t) in portfolio value $\eqref{portDiff2} expanded by the given stock price of the stochastic differential equation \eqref{stockPrice} can be interpreted as follows.

  • $rX(t)dt$: A term that reflects the profit from the expected rate of return of the portfolio $r$
  • $\Delta(t)( \alpha -r )S(t)dt$: A term that reflects the risk premium on equity investments $\alpha-r$
  • $\Delta(t)\sigma S(t)dW(t)$: Volatility term proportional to the size of equity investment

Then, consider the discounted stock price $\mathrm{e}^{-rt} S(t)$ and the discounted portfolio value $\mathrm{e}^{-rt}X(t)$. In order to analyze these price processes, it should be noted that they are different from ordinary analysis due to the nature of the geometric Brownian motion that appeared as a model of stock prices. This property means that Brownian motion has a non-zero quadratic variation, which, as it will be seen later, is the source of the volatility term of the Black-Scholes PDE. In other words, to differentiate an expression expressed in the form of $f(W(t))$ with $f(x)$ as a differentiable function and $W(t)$ as the Brown motion. If Brown motion $W(t)$ is differentiable, the following equation holds from ordinary analysis by the chain rule of the composite function.

\begin{eqnarray} \frac{d}{dt}f(W(t))=f’(W(t))W’(t) \Leftrightarrow df(W(t))=f’(W(t))W’(t)dt=f’(W(t))dW(t)． \end{eqnarray}

However, since Brownian motion accumulates a quadratic variation of 1 per unit time $(dWdW = dt)$, the derivative of $f(W(t))$ has a correction term shown in the following equation.

\begin{eqnarray}\label{ItoFormDiff} df(W(t))=f’(W(t))dW(t)+\frac{1}{2}f’‘(W(t))dt． \end{eqnarray}

The formula \eqref{ItoFormDiff} is called Ito’s formula in differential form. The equation \eqref{ItoFormDiff} is an equation obtained by differentiating Ito’s formula \eqref{ItoForm} for the Brownian motion shown below.

  • Ito’s formula for Brownian motion

Partial differentiation of $f(t, x)$, $f_ {t}(t, x)$, $f_{x}(t, x)$, $f_ {xx}(t, x)$ is defined and continuous. Let $W(t)$ be the Brownian motion. At this time, the following equation holds for all $T \geq0$.

\begin{eqnarray}\label{ItoForm} f(T，W(T))=f(0，W(0))+\int_{0}^{T}f_{t}(t，W(t))dt+\int_{0}^{T}f_{x}(t，W(t))dW(t)+\frac{1}{2}\int_{0}^{T}f_{xx}(t，W(t))dt． \end{eqnarray}

Where the third term of the expression \eqref{ItoForm} $\int_{0}^{T}f_{x}(t, W(t)) dW(t)$ is called Ito integral, and the second and fourth terms are Lebesgue integrals for variables of time. Ito’s formula can be extended to general stochastic processes including the Brown process, excluding jump process. This stochastic process is called the Ito process $dX(t)=a(t)dW(t)+b(t)dt$, which accumulate the quadratic variation of $a^{2}(t)$ per unit time. It is assumed that $a(t)$ and $b(t)$ are adapted stochastic processes.

  • Ito’s formula for the Ito process

Let $X(t)$ be the Ito process, and $f(t, x)$ be the partial derivative $f_{t}(t, x)$, $f_{x}(t, x)$, $f_{xx}(t, x)$ be a defined and continuous function, and let $W(t)$ be Brownian motion. At this time, the following equation holds for all $T \geq 0$.

\begin{eqnarray}\label{ItoFormII} f(T，X(T))=f(0，X(0))+\int_{0}^{T}f_{t}(t，X(t))dt+\int_{0}^{T}f_{x}(t，X(t))dX(t)+\frac{1}{2}\int_{0}^{T}f_{xx}(t，X(t))a^{2}(t)dt． \end{eqnarray}

Therefore, the derivative of the discounted stock price $\mathrm{e}^{-rt} S(t)$ is $f(t, x) = \mathrm{e}^{-rt} x$ which is gotten to apply Ito’s formula \eqref{ItoForm}. In the following, it is used that $f_{t} (t, x )=-r \mathrm{e}^{-rt} x$, $ f_{x} (t, x )=\, mathrm{e}^{-rt}$, $f_{xx} (t, x )=0 $.

\begin{eqnarray}\label{DiscountStockPrice} d( \mathrm{e}^{-rt}S(t) )=df(t，S(t)) =f_{t}(t，S(t))dt+f_{x}(t，S(t))dS(t)+f_{xx}(t，S(t))dt=-r\mathrm{e}^{-rt}S(t)dt+\mathrm{e}^{-rt}dS(t) =-r\mathrm{e}^{-rt}S(t)dt+\mathrm{e}^{-rt}( \alpha S(t)dt+\sigma S(t)dW(t))=(\alpha-r)\mathrm{e}^{-rt}S(t)dt+\sigma\mathrm{e}^{-rt}S(t)dW(t) \end{eqnarray}

Similarly, the derivative of the discounted portfolio value $\mathrm{e}^{-rt} X(t)$ is given by

\begin{eqnarray}\label{DiscountPort} d( \mathrm{e}^{-rt}X(t) )=df(t，X(t))=f_{t}(t，X(t))dt+f_{x}(t，X(t))dX(t)+\frac{1}{2}f_{xx}(t，X(t))dt=-r\mathrm{e}^{-rt}X(t)dt+\mathrm{e}^{-rt}dX(t) =-r\mathrm{e}^{-rt}X(t)dt+\mathrm{e}^{-rt}( \Delta(t)dS(t)+r(X(t)-\Delta(t)S(t))dt)=\Delta(t)(\alpha-r)\mathrm{e}^{-rt}S(t)dt+\Delta(t)\sigma\mathrm{e}^{-rt}S(t)dW(t)=\Delta(t)d(\mathrm{e}^{-rt}S(t))． \end{eqnarray}

In this way, the formula \eqref{DiscountPort} is obtained by the discounted stock price $S(t)$ and the portfolio value $X(t)$ at the discount rate $\mathrm{e}^{-rt}$ with the interest rate $r$ of the money market account. The formula \eqref{DiscountPort} indicates that small change in the discount portfolio value $d(\mathrm{e}^{-rt}X(t))$ depends only on a small change in the discounted stock price $d(\mathrm{e}^{-rt}S(t))$.

Consider a European call option where the payment on the maturity date $T$ is $(S(T)-K )^{+}$. Here, the strike price $K$ is assumed to be a non-negative constant. Fischer Black, Myron Scholes and Robert Cox Merton argued that the value of European call option on maturity $T$ with payment $(S(T) -K )^{+}$ at any time $t$ depends on the length of time to maturity, the value of the stock price $S(t)$ at time $t$, the model parameters $r, \sigma$ and the strike price $K$.

Since the model parameters $r, \sigma$ and the strike price $K$ are assumed to be constants, among the variables that affect the price evaluation of the European call option, only the time $t$ and the stock price $S(t)$ are changed.

Here, according to the inference of Fischer Black et al, let the value of the European call option be $c(t, x)$ when the stock price is $S(t) = x$ at time $t$. Since the stock price $S(t)$ follows the stochastic differential equation \eqref{stockPrice}, the value of the European call option is also stochastic.

Therefore, the value of the European call option is the stochastic process $c(t, S(t))$ obtained by replacing the stochastic process $S(t)$ with the parameter $x$. At the initial time, the future stock price $S(t)$ is unknown because it is probabilistic, then the value $c(t, S(t))$ of the future European call option is also unknown. Therefore, if the function $c(t, x)$ is determined, the value of the future European call option $c(t, S(t))$ is expressed by the future stock price $S(t)$.

To duplicate a portfolio which does full hedge the option sell position, it is required to fit the value of the duplicate portfolio $X(t)$ at any time $t\in[0, T]$ in the value of future European call options $c(t, S(t))$ starting with a certain initial fund of $X(0)$ and investing in the stock market and money market accounts.

The necessary and sufficient conditions for constructing this replication portfolio are $\ mathrm{e}^{-rt}X(t)=\mathrm{e}^{-rt} c(t, S (t))$ for any time $t \in[0, T]$. To show that the necessary and sufficient conditions are met,

\begin{eqnarray}\label{IsCopyPort} d(\mathrm{e}^{-rt}X(t))=d(\mathrm{e}^{-rt}c(t，S(t)))，t\in[0，T) \end{eqnarray}

And show that $X(0) = c(0, S(0))$ holds. In fact, the following equation is obtained by integrating the equation \eqref{IsCopyPort} from 0 to $t$.

\begin{eqnarray} \mathrm{e}^{-rt}X(t)-X(0)=\mathrm{e}^{-rt}c(t，S(t))-c(0，S(0))，t\in[0，T) \end{eqnarray}

At this time, if $X(0)=c(0,S(0))$, then it is obtained the necessary and sufficient conditional expression $\mathrm{e}^{-rt}X(t)=\mathrm{e}^{-rt}c(t,S(t))$.

Necessary and sufficient condition to duplicate the portfolio, which does fully hedge the option sell position

To expand the right-hand side \eqref{IsCopyPort}, it is required to calculate the derivative of discount option price $\mathrm{e}^{-rt} c(t, S (t, S) t))$. To do this, it is needed to calculate the derivative of the option price $c(t, S(t))$. From Ito’s formula for the Ito process \eqref{ItoFormII} and the given stock price stochastic differential equation \eqref{stockPrice}

\begin{eqnarray} dc(t，S(t))=c_{t}(t，S(t))dt+c_{x}(t，S(t))dS(t)+\frac{1}{2}c_{xx}(t，S(t))(\sigma S(t))^{2}dt=c_{t}(t，S(t))dt+c_{x}(t，S(t))( \alpha S(t)dt+\sigma S(t)dW(t) )+\frac{1}{2}c_{xx}(t，S(t))(\sigma S(t))^{2}dt=( c_{t}(t，S(t))dt+\alpha c_{x}(t，S(t))+\frac{1}{2}\sigma^{2} S^{2}(t)c_{xx}(t，S(t)))dt+\sigma S(t)c_{x}(t，S(t))dW(t)． \end{eqnarray}

It is calculated the derivative of discount option price $\mathrm{e}^{-rt} c (t)$ by Ito’s formula \eqref{ItoFormII} for the Ito process as $f(t, x)=\mathrm{e}^{-rt}x,S(t))$. In the following, it is used that $f_{t}(t,x)=-r\mathrm{e}^{-rt}x$, $f_{x}(t,x)=\mathrm{e}^{-rt}$, $f_{xx}(t,x)=0$.

\begin{eqnarray}\label{expandCallDis} d(\mathrm{e}^{-rt}c(t，S(t)))=df(t，c(t，S(t)))=f_{t}(t，c(t，S(t)))dt+f_{x}(t，c(t，S(t)))dc(t，S(t))+\frac{1}{2}f_{xx}(t，c(t，S(t)) (\sigma S(t)c_{x}(t，S(t)))^{2}dt=-r\mathrm{e}^{-rt}c(t，S(t))+\mathrm{e}^{-rt} dc(t，S(t))=\mathrm{e}^{-rt}(-rc(t，S(t))+c_{t}(t，S(t))+\alpha c_{x}(t，S(t))+\frac{1}{2}\sigma^{2} S^{2}(t)c_{xx}(t，S(t)))dt+\mathrm{e}^{-rt}\sigma S(t)c_{x}(t，S(t))dW(t)． \end{eqnarray}

It is given the following equation to expand the right-hand side of \eqref{IsCopyPort} with the formula \eqref{expandCallDis} and the left-hand side with the formula \eqref{DiscountPort}.

\begin{eqnarray}\label{expandCopyPort} \Delta(t)(\alpha-r)S(t)dt+\Delta(t)\sigma S(t)dW(t)=( -rc(t，S(t))+c_{t}(t，S(t))+\alpha c_{x}(t，S(t))+\frac{1}{2}\sigma^{2} S^{2}(t)c_{xx}(t，S(t)))dt+\sigma S(t)c_{x}(t，S(t))dW(t)． \end{eqnarray}

If the expansion formula \eqref{expandCopyPort}, which is a necessary and sufficient condition for replicating a portfolio that completely hedges an option’s sell position, holds, the coefficients on both sides must be matched. Therefore, since the coefficients on both sides of $dW(t)$ are equal, the following equation is obtained.

\begin{eqnarray}\label{deltaHedging} \Delta(t)=c_{x}(t，S(t))，t\in[0，T) \end{eqnarray}

The number of shares held to fully hedge the option’s sell position at any time $t$ before maturity is the partial derivative of the option value at time $t$ with respect to the stock price $c_{x}(t, S(t))$. The relational expression \eqref{deltaHedging} of $\Delta(t)$ is called the delta-hedging rule because it is called the optional delta.

Since the coefficients on both sides of $dt$ are also equal, the following equation holds for any time $t\in[0, T)$ from the relational expression \eqref{deltaHedging} of $\Delta(t)$.

\begin{eqnarray}\label{expandCopyPortForT} (\alpha-r)S(t)c_{x}(t，S(t)) =-rc(t，S(t))+c_{t}(t，S(t))+\alpha c_{x}(t，S(t))+\frac{1}{2}\sigma^{2} S^{2}(t)c_{xx}(t，S(t)) \end{eqnarray}

Since $\alpha c_{x}(t,S(t))$ term appears on both sides of the expression \eqref{expandCopyPortForT} for any time $t\in[0, T)$, it is obtained the following equation to extract this term.

\begin{eqnarray}\label{PreBlackScholes} rc(t，S(t))=c_{t}(t，S(t))+rS(t)c_{x}(t，S(t))+\frac{1}{2}\sigma^{2} S^{2}(t)c_{xx}(t，S(t)) \end{eqnarray}

Here, if the value of the European call option when the stock price is $S(t) = x$ at time $t$ is $c(t, x)$, then it is obtained the partial differential equation of Black-Scholes from the equation \eqref{PreBlackScholes}.

\begin{eqnarray}\label{BlackScholes} c_{t}(t，x)+rxc_{x}(t，x)\frac{1}{2}\sigma^{2} x^{2}c_{xx}(t，x)=rc(t，x)，t\in[0，T)，x\geq0 \end{eqnarray}

In addition to satisfying the equation \eqref{BlackScholes}, the function $c(t, x)$ must also satisfy the boundary conditions shown in the following equation.

\begin{eqnarray}\label{BoundaryCondtion} c(T，x)=(x-K)^{+}，c(t，0)=0，\lim_{x \to \infty}[c(t，x)-(x-\mathrm{e}^{-r( T-t )}K)]=0 \end{eqnarray}

In fact, the first equation of the expression \eqref{BoundaryCondtion} follows from the definition of the European call option.

The second equation is that when $x=0$ is substituted into the Black-Scholes partial differential equation \eqref{BlackScholes}, it is obtained the ordinary differential equation $c_{t}(t,0)=rc(t,0)$ for the function $c(t, 0)$ at time $t$ and then it is obtained its solution $c(t,0)=\mathrm{e}^{rt}$ is substituted $t=T$ for $c(0,0)$, which results in $c(T,0)=(0-K)^{+}=0$, $C(0,0)=0$.

In the third equation, ‘Put/Call Parity’ holds because it is likely to end with ‘In the money’ for European call options for very large $x$ and the value of the forward contract is almost equal to $x-\mathrm{e}^{-r (Tt )}K$.

In order to obtain an analytical solution of the Black-Scholes PDE \eqref{BlackScholes}, it is reduced by change of variables to a PDE whose exact solution is known. Therefore, two variables $u$ and $v$ are defined by the following equation, which are introduced for this purpose.

\begin{eqnarray}\label{val1} u=\log \frac{x}{K}+(r-\frac{\sigma^{2}}{2})(T-t) v=T-t \end{eqnarray}

Here, using the variables $u$ and $v$ defined in the expression \eqref{val1}, the function $c(t, x)$ is converted into two functions $\mathrm{e}^{-rv}$ and $y (u, v )$ as shown in the following expression.

\begin{eqnarray}\label{val2} c(t，x)=\mathrm{e}^{-rv}y(u，v)． \end{eqnarray}

Substituting the equation \eqref{val2} into the partial differential equation \eqref{BlackScholes} of Black-Scholes and rearranging it gives the following equation.

\begin{eqnarray}\label{HeatConductionEquation} y_{u u}(u，v)-\frac{2}{\sigma^{2}}y_{v}(u，v)=0． \end{eqnarray}

The equation \eqref{HeatConductionEquation} is a known partial differential equation with an exact solution called the heat conduction equation. Here, the boundary condition is changed by variable transformation.

$c(t，x)=$ $x-K\cdots x\geq K$ $0\cdots x \lt K$

Note that the following conditions have changed.

$y(u，0)=$ $K( \mathrm{e}^{u}-1 )\cdots u\geq 0$ $0\cdots u \lt 0$

It is known that the heat conduction equation \eqref{HeatConductionEquation} is solved by the separation of variables method. The separation of variables method means that the solution $y(u,v)$ is expressed as the product of two functions $p(u)$ and $q(v)$.

\begin{eqnarray}\label{val3} y(u，v)=p(u)q(v)． \end{eqnarray}

Substituting the equation \eqref{val3} into the heat conduction equation \eqref{HeatConductionEquation} and rearranging it gives the following equation.

\begin{eqnarray}\label{SeparationVariables} \frac{p_{u u}(u)}{p(u)}=\frac{2}{\sigma^{2}}\frac{q_{v}(v)}{q(v)}． \end{eqnarray}

Focusing on both sides of the obtained equation \eqref{SeparationVariables}, the left side is a function of only $u$ and the right side is a function of only $v$, then constants irrelevant to $u$ and $v$, which is $-k^{2}(0\leq k\lt \infty)$. Therefore, finding the solution of the heat conduction equation \eqref{HeatConductionEquation} results in the problem of solving the following two ordinary differential equations.

  • $\dfrac{p_{u u}(u)}{p(u)}=-k^{2}$
  • $\dfrac{2}{\sigma^{2}}\dfrac{q_{v}(v)}{q(v)}=-k^{2}$

The first equation is a second-order ordinary differential equation with constant coefficients, then if the coefficients are constants $C(k)$ and $D(k)$ with respect to constants $k$ irrelevant to $u$ and $v$, the solution is given by the following equation.

\begin{eqnarray}\label{ODE1Sol} p(u)=C(k)\sin( ku ) +D(k)\cos( ku ) ． \end{eqnarray}

Since the second equation is a variable separable ordinary differential equation, if the coefficients are the constant $E(k)$ for the constant $k$ irrelevant to $u$ and $v$, the solution is given by the following equation.

\begin{eqnarray}\label{ODE2Sol} q(v)=E(k)\exp( -\frac{\sigma^{2}k^{2}}{2}v ) ． \end{eqnarray}

Therefore, the solution $y (u, v )$ of the heat conduction equation \eqref{HeatConductionEquation} is a constant $k (0\leq k \lt \infty )$ irrelevant to $u$ and $v$. Then all are superposed to obtain the general solution shown in the following equation.

\begin{eqnarray}\label{HeatConductionEquationSol} y(u，v)=\int_{0}^{+\infty}( C(k)\sin( ku ) +D(k)\cos( ku ) )\exp( -\frac{\sigma^{2}k^{2}}{2}v )dk． \end{eqnarray}

Here, the coefficients $C(k)$ and $D(k)$ of the equation \eqref{HeatConductionEquationSol} are given by the following equation from Fourier’s integral theorem as $y(u,0)=g(a)$.

\begin{eqnarray}\label{FourierCoefficient} C(k)=\frac{1}{\pi}\int_{-\infty}^{+\infty}g(a)\cos (ka) da \end{eqnarray} \begin{eqnarray}\label{FourierCoefficient2} D(k)=\frac{1}{\pi}\int_{-\infty}^{+\infty}g(a)\sin (ka) da \end{eqnarray}

Substituting the obtained equation of coefficients \eqref{FourierCoefficient} and \eqref{FourierCoefficient2} into the solution \eqref{HeatConductionEquationSol} of the heat conduction equation \eqref{HeatConductionEquation} and rearranging it gives the following equation.

\begin{eqnarray}\label{HeatConductionEquationSol2} y(u，v)=\dfrac{1}{\sigma \sqrt[]{\mathstrut 2 \pi v}} \int_{-\infty}^{+\infty}g(a)\exp( -\dfrac{1}{2}(\dfrac{a-u}{\sigma \sqrt[]{\mathstrut v}})^{2} ) da． \end{eqnarray}

For further analysis, the change of variables is performed again as follows. \begin{eqnarray}\label{} w=\dfrac{a-u}{\sigma \sqrt[]{\mathstrut v}}． \end{eqnarray}

Note that the boundary conditions have changed to the following conditions due to the change of variables.

$g(a)=g( u+\sigma \sqrt[]{\mathstrut v}w )=$

$K( \mathrm{e}^{u+\sigma \sqrt[]{\mathstrut v}w}-1 )\cdots w\geq -\dfrac{u}{\sigma \sqrt[]{\mathstrut v}}$

Or $0\cdots w \lt -\dfrac{u}{\sigma \sqrt[]{\mathstrut v}}$

Focusing on the change of variables by $w$ and the changed boundary conditions, the equation \eqref{HeatConductionEquationSol2} is expanded as follows.

\begin{eqnarray}\label{HeatConductionEquationSol3} y(u，v)=\dfrac{1}{\sqrt[]{\mathstrut 2 \pi}} \int_{-\dfrac{u}{\sigma \sqrt[]{\mathstrut v}}}^{+\infty}K( \exp(u+\sigma \sqrt[]{\mathstrut v}w)-1 )\exp( -\dfrac{w^{2}}{2} )dw =\dfrac{1}{\sqrt[]{\mathstrut 2 \pi}} \int_{-\frac{u}{\sigma \sqrt[]{\mathstrut v}}}^{+\infty}K\exp(u+\sigma \sqrt[]{\mathstrut v}w-\dfrac{w^{2}}{2}) dw-\dfrac{1}{\sqrt[]{\mathstrut 2 \pi}} \int_{-\frac{u}{\sigma \sqrt[]{\mathstrut v}}}^{+\infty} K \exp(-\dfrac{w^{2}}{2})dw． \end{eqnarray}

Furthermore, if the first term on the right side of the equation \eqref{HeatConductionEquationSol3} is replaced with $z = w- \sigma \sqrt[]{\mathstrut v}$ and the second term is replaced with $z = w$, the following equation is obtained.

\begin{eqnarray}\label{HeatConductionEquationSol4} y(u，v)=x\exp(rv) \dfrac{1}{ \sqrt[]{\mathstrut 2 \pi}} \int_{-\frac{u}{\sigma \sqrt[]{\mathstrut v}}^{+\infty}-\sigma \sqrt[]{\mathstrut v}}\exp( -\dfrac{z^{2}}{2} ) dz-K\dfrac{1}{ \sqrt[]{\mathstrut 2 \pi}} \int_{-\frac{u}{\sigma \sqrt[]{\mathstrut v}}}^{+\infty}\exp( -\dfrac{z^{2}}{2} ) dz =x\exp(rv)N(\frac{u}{\sigma \sqrt[]{\mathstrut v}}+\sigma \sqrt[]{\mathstrut v})-KN( \frac{u}{\sigma \sqrt[]{\mathstrut v}})． \end{eqnarray}

Finally, by substituting the obtained equation \eqref{HeatConductionEquationSol4} into the solution of the Black-Scholes PDE \eqref{val2} is obtained as an analytical solution of the Black-Scholes PDE, which is a function that satisfies the Black-Scholes PDE.

\begin{eqnarray}\label c(t，x)=xN(\frac{u}{\sigma \sqrt[]{\mathstrut v}} + \sigma \sqrt[]{\mathstrut v} ) -K \exp(-rv)N(\frac{u}{\sigma \sqrt[]{\mathstrut v}}) \end{eqnarray}

However,

$u=\log \frac{x}{K}+(r-\frac{\sigma^{2}}{2})(T-t)$

$v=T-t$

Reference

  • Shreve, S., 2004. Stochastic Calculus for Finance II: Continuous-Time Models
  • Ishimura, S., 2008. Black Shoals differential equations for finance and securities (In Japanese)