{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![MOSEK ApS](https://www.mosek.com/static/images/branding/webgraphmoseklogocolor.png )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Utility based option pricing with transaction costs and diversification"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this notebook, we present the Fusion implementation of the utility based option pricing model, presented by Andersen et. al. (1999) (https://www.sciencedirect.com/science/article/pii/S0168927498001044). The purpose of the model is to estimate the reservation purchase price of a European call option, written on a risky security when there is proportional transaction costs in the market. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The Model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider an economy evolving over T periods, ${0,t_1,t_2,t_3,...,t_T = \\overline{T}}$ where $\\overline{T}$ is the time horizon in years. We consider that the investor has one risk-free security, such as a bond, that pays at a constant interest rate of $r\\geq 0$, and $m$ risky securities with price processes denoted by $(S_1,S_2,...,S_m)$. The price of the risky securities evolves in an event tree such that: \n",
"\n",
"$$(S_{1,n},S_{2,n},..,S_{m,n}) = \\begin{cases}\n",
"(a_1 S_{1,n-1}, a_2 S_{2,n-1},..,a_m S_{m,n-1}), \\,\\, \\text{with probability } q_1, \\\\\n",
"(b_1 S_{1,n-1}, b_2 S_{2,n-1},..,b_m S_{m,n-1}), \\,\\, \\text{with probability } q_2, \\\\\n",
"(c_1 S_{1,n-1}, c_2 S_{2,n-1},..,c_m S_{m,n-1}), \\,\\, \\text{with probability } q_3 = (1-q_1-q_2)\n",
"\\end{cases}\n",
"$$\n",
"\n",
"\n",
"where the three possibilities lead to an event tree of splitting index 3 (In the fusion model presented below, we have considered a general case where one can have a different splitting index). The event tree for two risky securities can be visualized as shown in the figure below."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let, $\\mathbf{\\alpha_{n}} = (\\alpha_1,\\alpha_2,...,\\alpha_m)$ denote the number of units of the risky security held by the investor at time $t_n$, and $\\beta_n$ denote the dollar amount held by the investor in bonds at the same time. Then, the 'budget constraint' on the bonds an investor can buy in the next period is given by:\n",
"\n",
"\n",
"$$\\beta_n(I_n) \\leq (1+r)\\beta_{n-1}(I_n) + \\sum_{i=1}^{m} S_{i,n}(I_n)[\\alpha_{i,n-1}(I_n) - \\alpha_{i,n}(I_n) - \\theta_i |\\alpha_{i,n-1}(I_n) - \\alpha_{i,n}(I_n)|]$$\n",
"\n",
"\n",
"where $\\theta_i$ is the transaction cost for trading of risky security $i$ and $I_n$ denotes the path being considered (in other words, the sequence of events in the event tree). The total number of paths possible for a tree of splitting index $s$, over $n$ time periods is $s^n$. Initially, if we consider that the wealth of the investor is $W_0$, then of course the first budget constraint becomes:\n",
"\n",
"$$\\beta_0 \\leq W_0 - \\sum_{i=1}^{m} S_{i,0}[\\alpha_{i,0} + |\\alpha_{i,0}|]$$\n",
"\n",
"\n",
"Additionally, in the final period the investor will sell all of the risky securities, thus ending up with a wealth $W_T(I_T)$ for path $I_T$, such that:\n",
"\n",
"$$W_T(I_T) \\leq (1+r)\\beta_{T-1}(I_T) + \\sum_{i=1}^{m} S_{i,T}(I_T)[\\alpha_{i,T-1}(I_T) - \\theta_i |\\alpha_{i,T-1}(I_T)|]$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 1.) Portfolio choice"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The goal of the investor is to choose a portfolio strategy (i.e. the sequence {${\\mathbf{\\alpha_n(I_n)},\\beta_n(I_n)}$}) that maximizes the expected utility. The expected utility is given by:\n",
"\n",
"$$E[U(W_T)] = \\sum_{I_{T}\\in F_T} q_1^{A(I_T)}q_2^{B(I_T)}(1-q_1-q_2)^{T-A(I_T)-B(I_T)}U(W_T(I_T))$$\n",
"\n",
"\n",
"Note that the summation is over all the possible paths for a tree of a given split index and T time-periods (the set denoted by $F_T$). $A(I_T)$ and $B(I_T)$ denote the number of times the first and second possibilities are considered in every path, respectively. Following the equations presented above, the complete optimization problem can be written as:\n",
"\n",
"\n",
"$$U^*(W_0) = \\text{maximize}_{({\\mathbf{\\alpha}_n(I_n),\\beta_n(I_n))}_{I_n\\in F_T,{n=1,2,..,T-1}}}\\,\\,\\,E[U(W_T)]$$\n",
"\n",
"\n",
"$$\\text{subject to: } \\beta_0 \\leq W_0 - \\sum_{i=1}^{m} S_{i,0}[\\alpha_{i,0} + |\\alpha_{i,0}|],$$\n",
"$$\\beta_n(I_n) \\leq (1+r)\\beta_{n-1}(I_n) + \\sum_{i=1}^{m} S_{i,n}(I_n)[\\alpha_{i,n-1}(I_n) - \\alpha_{i,n}(I_n) - \\theta_i |\\alpha_{i,n-1}(I_n) - \\alpha_{i,n}(I_n)|],$$\n",
"$$W_T(I_T) \\leq (1+r)\\beta_{T-1}(I_T) + \\sum_{i=1}^{m} S_{i,T}(I_T)[\\alpha_{i,T-1}(I_T) - \\theta_i |\\alpha_{i,T-1}(I_T)|],\\,\\,\\, \\forall I_T\\in F_T$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 2.) Price vector process"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The continuous price process that we intended to approximate is:\n",
"\n",
"$$\\frac{d\\tilde{S_i}}{\\tilde{S_i}} = \\mu_idt + \\sum_{j=1}^{m}\\sigma_{ij} dW_j \\,\\, , \\,\\, \\forall i = 1,2,...,m$$\n",
"\n",
"where $\\mu_i$ and $\\sigma_{ij}$ are positive constants and the $W_j$ denote un-correlated standard Wiener processes. To construct a discrete approximation, we first define a stochastic vector as follows:\n",
"\n",
"$$(\\epsilon_1,\\epsilon_2,..,\\epsilon_m) = \\begin{cases}\n",
"(\\epsilon_1 (\\omega_1),\\epsilon_2 (\\omega_1),..,\\epsilon_m (\\omega_1)), \\,\\, \\text{with probability } q_1, \\\\\n",
"(\\epsilon_1 (\\omega_2),\\epsilon_2 (\\omega_2),..,\\epsilon_m (\\omega_2)), \\,\\, \\text{with probability } q_2, \\\\\n",
"(\\epsilon_1 (\\omega_3),\\epsilon_2 (\\omega_3),..,\\epsilon_m (\\omega_3)), \\,\\, \\text{with probability } q_3 = (1-q_1-q_2)\n",
"\\end{cases}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The equation describing the event tree (first equation of the notebook) becomes a discrete approximation of the above-stated stochastic process if we set:\n",
"\n",
"$$a_i = e^{\\mu_i \\Delta t}\\bigg( \\frac{e^{\\big[\\sum_{j=1}^{m}\\sigma_{ij}\\epsilon_j(\\omega_1)\\sqrt{\\Delta t}\\big]}}{z_i}\\bigg)$$\n",
"\n",
"$$b_i = e^{\\mu_i \\Delta t}\\bigg( \\frac{e^{\\big[\\sum_{j=1}^{m}\\sigma_{ij}\\epsilon_j(\\omega_2)\\sqrt{\\Delta t}\\big]}}{z_i}\\bigg)$$\n",
"\n",
"$$c_i = e^{\\mu_i \\Delta t}\\bigg( \\frac{e^{\\big[\\sum_{j=1}^{m}\\sigma_{ij}\\epsilon_j(\\omega_3)\\sqrt{\\Delta t}\\big]}}{z_i}\\bigg)$$\n",
"\n",
"\n",
"\n",
"In the limit $T\\to \\infty$ (or alternatively $\\Delta t \\to 0$), the approximation approaches the continuous process."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3.) Investor's reservation purchase price for a European call option"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Suppose that an investor is interested in buying a European call option on the security $S_1$, with time to maturity $\\overline{T}$ and a strike price $K>0$. At the expiration time, the investor would get a payment of $\\text{max}(S_{1,T} - K, 0)$ (assuming cash settlement and also that the investor will not be re-selling the option once it has been purchased.) Our goal is now to estimate the highest price this investor is willing to pay, to purchase such an option. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If the investor does not purchase the call option and only trades in the risky securities $S_i$ and the bonds, then the portfolio is given by the optimization problem stated above. However, if the investor buys the call option at a reservation purchase price of $C$, then their portfolio becomes the following:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$P(C,W_0) = \\text{maximize}_{({\\mathbf{\\alpha}_n(I_n),\\beta_n(I_n))}_{I_n\\in F_T,{n=1,2,..,T-1}}}\\,\\,\\,E[U(W_T)]$$\n",
"\n",
"\n",
"$$\\text{subject to: }\\beta_0 \\leq W_0 - C - \\sum_{i=1}^{m} S_{i,0}[\\alpha_{i,0} + |\\alpha_{i,0}|],$$\n",
"$$\\beta_n(I_n) \\leq (1+r)\\beta_{n-1}(I_n) + \\sum_{i=1}^{m} S_{i,n}(I_n)[\\alpha_{i,n-1}(I_n) - \\alpha_{i,n}(I_n) - \\theta_i |\\alpha_{i,n-1}(I_n) - \\alpha_{i,n}(I_n)|],$$\n",
"$$W_T(I_T) \\leq (1+r)\\beta_{T-1}(I_T) + \\sum_{i=1}^{m} S_{i,T}(I_T)[\\alpha_{i,T-1}(I_T) - \\theta_i |\\alpha_{i,T-1}(I_T)|] + \\text{max}(S_{1,T}(I_T) - K,0), \\,\\,\\, \\forall I_T\\in F_T$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The highest price that investor is willing to pay for the given call option, is therefore given by the price for which the maximum expected utility in the two portfolios becomes equal, making the investor indifferent to the choices. Thus, the final optimization problem that we need to consider is:\n",
"\n",
"$$C^* = \\text{maximize}_{({\\mathbf{\\alpha}_n(I_n),\\beta_n(I_n))}_{I_n\\in F_T,{n=1,2,..,T-1}}} \\,\\,\\, C$$\n",
"\n",
"$$\\text{subject to: } E[U(W_T)] \\geq U^*(W_0) $$\n",
"\n",
"along with the other constraints mentioned in the optimization problem for the portfolio in which the option is purchased."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 4.) Utility function"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the optimization problems mentioned above, we need to optimize the expected utility. The utility functions that we consider are members of a family called HARA utility functions (Hyperbolic Absolute Risk Aversion). The general expression for HARA utility is:\n",
"\n",
"$$U(W) = \\frac{1-\\gamma}{\\gamma}\\bigg(\\frac{aW}{1-\\gamma} + b\\bigg)^\\gamma \\,;\\, a>0$$\n",
"\n",
"where $W > ((\\gamma - 1)b)/a$. Note that the exponential and logarithmic utility functions are also members of the HARA-class."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Another thing to consider is the Arrow-Pratt measure of the Absolute risk aversion. For the HARA-class, it is:\n",
"\n",
"$$ARA(W) = \\bigg(\\frac{W}{1-\\gamma} + \\frac{b}{a}\\bigg)^{-1}$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fusion Implementation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we proceed with the construction of the fusion model. We start by making a Tree class that will represent the event tree discussed above."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The various constants in the initialization of the Tree class are as follows:\n",
" M : Model, T : Time steps, W_0 : Initial wealth, S : Price-scaling matrix, theta : Transaction costs, S_v_0 : Initial price vector for the risky securities, r : Interest rate for the bonds, U : [a,b,$\\gamma$] for the HARA utility parameters ."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### HARA utility as a Power cone:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The objective function is:\n",
"\n",
"$$\\text{maximize }E[U(W_T)] = \\sum_{I_T\\in F_T}P(I_T) \\bigg(\\frac{1-\\gamma}{\\gamma}\\bigg) \\bigg(\\frac{aW(I_T)}{1-\\gamma} + b\\bigg)^\\gamma$$\n",
"\n",
"where $P(I_T)$ is the probability of a given path. We can re-write this as:\n",
"\n",
"\n",
"$$\\text{maximize } \\sum_{I_T\\in F_T}P(I_T) \\bigg(\\frac{1-\\gamma}{\\gamma}\\bigg) h(I_T)$$\n",
"\n",
"$$\\text{subject to : } h(I_T)\\leq \\bigg(\\frac{aW(I_T)}{1-\\gamma} + b\\bigg)^\\gamma$$\n",
"\n",
"\n",
"Now, the constraint can be expressed as a Power cone. However, there are a few possible cases:\n",
"\n",
"1.) $\\gamma > 1$: In this case, the conic representation is:\n",
"\n",
"$$\\Bigg(h(I_T),1,\\bigg(\\frac{aW(I_T)}{1-\\gamma} + b\\bigg) \\Bigg) \\in \\mathcal{P}_3^{1/\\gamma}$$\n",
"\n",
"2.) $0 < \\gamma < 1$:\n",
"\n",
"$$\\Bigg(\\bigg(\\frac{aW(I_T)}{1-\\gamma} + b\\bigg),1, h(I_T) \\Bigg) \\in \\mathcal{P}_3^{\\gamma}$$\n",
"\n",
"3.) $\\gamma < 0$:\n",
"\n",
"$$\\Bigg(h(I_T),\\bigg(\\frac{aW(I_T)}{1-\\gamma} + b\\bigg),1 \\Bigg) \\in \\mathcal{P}_3^{1/(1-\\gamma)}$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exponential utility as an exponential cone:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The objective function in the case of exponential utility will be given by:\n",
"\n",
"$$\\text{maximize }E[U(W_T)] = \\sum_{I_T\\in F_T}P(I_T) \\bigg(1 - e^{- [\\text{ARA}(W)] W(I_T)}\\bigg)\n",
"$$\n",
"\n",
"if the Absolute Risk aversion is nonzero. However, if ARA$(W) = 0$, then the objective function is simply:\n",
"\n",
"$$\\text{maximize }E[U(W_T)] = \\sum_{I_T\\in F_T}P(I_T) W(I_T)$$\n",
"\n",
"For the non-zero ARA$(W)$, we express the objective function as follows:\n",
"\n",
"$$\\text{maximize } \\sum_{I_T\\in F_T}P(I_T) h(I_T)$$\n",
"\n",
"$$\\text{subject to: } (1-h(I_T))\\geq e^{-[\\text{ARA}(W)]W(I_T)}$$\n",
"\n",
"and the constraint is readily expressed as an exponential cone:\n",
"\n",
"$$\\big((1-h(I_T)),1,-[\\text{ARA}(W)]W(I_T)\\big) \\in K_{exp}$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from mosek.fusion import *\n",
"import sys\n",
"import matplotlib.pyplot as plt\n",
"import time"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"class Tree: \n",
" #Instantiate the Tree class.\n",
" def __init__(self,M,T,W_0,S,theta,S_v_0,r,U):\n",
" self.M = M\n",
" self.W_0 = W_0\n",
" self.T = T\n",
" self.S_v = S\n",
" self.r = r \n",
" #Shape of the cost scaling matrix will give us the split index and the number of risky securities.\n",
" self.s,self.n = S.shape\n",
" self.cost = np.asarray([theta])\n",
" self.S_v_0 = S_v_0\n",
" #Parameters for the utility function.\n",
" self.U = U\n",
" self.S_v_t = S_v_0\n",
" #This is to decide which utility will be used.\n",
" self.util_dispatch = {'HARA': self.HARA_util,'EXP': self.exp_util}\n",
" self.levels = []\n",
" \n",
" \n",
" #Method to make all the levels, or the time steps in the tree.\n",
" def level_make(self):\n",
" #Number of risky securities, i.e. alpha{i,t=0} [Note that self.n is the number of risky securities]\n",
" self.a0 = self.M.variable('a_0',[1,self.n])\n",
" #Bonds at t=0.\n",
" self.b0 = self.M.variable('b_0',1)\n",
" #Making the first budget constraint.\n",
" self.budg_ex = Expr.sub(Expr.sub(self.W_0,self.b0),\n",
" Expr.dot(self.S_v_0,Expr.mulElm(1+self.cost,self.a0)))\n",
" self.root_budget = self.M.constraint('Budget_0',self.budg_ex,Domain.greaterThan(0.0))\n",
" #This list will hold the level objects associated to this tree. (Please look at the sub-class: Level)\n",
" a,b = self.a0,self.b0\n",
" for i in range(1,self.T+1):\n",
" #Appending level n, based on the level (n-1).\n",
" self.levels.append(Level(i,self,a,b))\n",
" a = self.levels[i-1].a_new\n",
" b = self.levels[i-1].b_new\n",
" #a_T is zero (). \n",
"\n",
" \n",
" #Method to make the HARA utility function and the corresponding Power-cone constraint.\n",
" def HARA_util(self,W):\n",
" self.h = self.M.variable('h',self.s**self.T,Domain.greaterThan(0.0))\n",
" \n",
" #HARA utility is only defined for when W > (gamma - 1)a/b.\n",
" self.M.constraint(W,Domain.greaterThan(self.U[1]*(self.U[2]-1)/self.U[0]))\n",
"\n",
" I = Expr.constTerm([self.s**self.T],1.0)\n",
" self.E1 = Expr.add(Expr.mul(self.U[0]/(1-self.U[2]),W),Expr.constTerm(W.getShape(),self.U[1]))\n",
" \n",
" #There are multiple cases for different values of gamma, modelled as follows\n",
" #For gamma > 1:\n",
" if self.U[2]>1.0:\n",
" self.M.constraint('Obj_terms_HARA',Expr.hstack(self.h,I,self.E1),Domain.inPPowerCone(1/self.U[2]))\n",
" else:\n",
" if self.U[2]<0.0: #For gamma < 0\n",
" self.M.constraint('Obj_terms_HARA',Expr.hstack(self.h,self.E1,I),Domain.inPPowerCone(1/(1-self.U[2])))\n",
" else: #For 0 < gamma < 1\n",
" self.M.constraint('Obj_terms_HARA',Expr.hstack(self.E1,I,self.h),Domain.inPPowerCone(self.U[2])) \n",
" return(Expr.mul(((1-self.U[2])/self.U[2]),self.h))\n",
" \n",
" \n",
" #Method to make the Exponential utility constraints.\n",
" def exp_util(self,W):\n",
" #Zero ARA(W) case\n",
" if self.U == 0:\n",
" return(W)\n",
" #Nonzero ARA(W) case\n",
" else:\n",
" self.h = self.M.variable('h',self.s**self.T)\n",
" I = Expr.constTerm([self.s**self.T],1.0)\n",
" E_exp = Expr.mul(-self.U,W)\n",
" self.M.constraint('Obj_terms_exp',Expr.hstack(Expr.sub(I,Expr.mul(self.U,self.h)),I,E_exp),\n",
" Domain.inPExpCone())\n",
" return(self.h)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Tree class defined above has the sub-class Level. The budget constraints involved in our model connect variables that lie in consecutive levels. Therefore, we can visualize all the constraints between two levels in the following manner (Note: for ease of visualization, we take the split index as $s = 3$, and a time period of n = 3):\n",
"\n",
"$$\n",
"\\begin{pmatrix}\n",
"a\\begin{bmatrix}\n",
"S_{(1,2)} \\\\\n",
"S_{(2,2)} \\\\\n",
"S_{(3,2)} \\\\\n",
"\\end{bmatrix}\\\\\n",
"b\\begin{bmatrix}\n",
"S_{(1,2)} \\\\\n",
"S_{(2,2)} \\\\\n",
"S_{(3,2)} \\\\\n",
"\\end{bmatrix}\\\\\n",
"c\\begin{bmatrix}\n",
"S_{(1,2)} \\\\\n",
"S_{(2,2)} \\\\\n",
"S_{(3,2)} \\\\\n",
"\\end{bmatrix}\\\\\n",
"\\end{pmatrix}\n",
":\n",
"\\begin{pmatrix}\n",
"\\begin{bmatrix}\n",
"\\mathbf{\\alpha_{(1,2)}},\\beta_{(1,2)} \\\\\n",
"\\mathbf{\\alpha_{(2,2)}},\\beta_{(2,2)} \\\\\n",
"\\mathbf{\\alpha_{(3,2)}},\\beta_{(3,2)} \\\\\n",
"\\end{bmatrix} \\\\\n",
"\\begin{bmatrix}\n",
"\\mathbf{\\alpha_{(1,2)}},\\beta_{(1,2)} \\\\\n",
"\\mathbf{\\alpha_{(2,2)}},\\beta_{(2,2)} \\\\\n",
"\\mathbf{\\alpha_{(3,2)}},\\beta_{(3,2)} \\\\\n",
"\\end{bmatrix} \\\\\n",
"\\begin{bmatrix}\n",
"\\mathbf{\\alpha_{(1,2)}},\\beta_{(1,2)} \\\\\n",
"\\mathbf{\\alpha_{(2,2)}},\\beta_{(2,2)} \\\\\n",
"\\mathbf{\\alpha_{(3,2)}},\\beta_{(3,2)} \\\\\n",
"\\end{bmatrix} \\\\\n",
"\\end{pmatrix} \\longrightarrow\n",
"\\begin{pmatrix}\n",
"\\mathbf{\\alpha_{(1,3)}},\\beta_{(1,3)} \\\\\n",
"\\mathbf{\\alpha_{(2,3)}},\\beta_{(2,3)} \\\\\n",
"\\mathbf{\\alpha_{(3,3)}},\\beta_{(3,3)} \\\\\n",
"\\mathbf{\\alpha_{(4,3)}},\\beta_{(4,3)} \\\\\n",
"\\mathbf{\\alpha_{(5,3)}},\\beta_{(5,3)} \\\\\n",
"\\mathbf{\\alpha_{(6,3)}},\\beta_{(6,3)} \\\\\n",
"\\mathbf{\\alpha_{(7,3)}},\\beta_{(7,3)} \\\\\n",
"\\mathbf{\\alpha_{(8,3)}},\\beta_{(8,3)} \\\\\n",
"\\mathbf{\\alpha_{(9,3)}},\\beta_{(9,3)} \\\\\n",
"\\end{pmatrix}\n",
"$$\n",
"\n",
"Here, the first column vector is a \"vstack\" of the price vector at level $n=2$, multiplied by the scaling factors for the price. Therefore, the first column is the new price vector (for $n=3$). In the second column, $\\alpha$'s will be vectors if there are multiple risky securities. Also, $\\alpha_{(i,n)}$ and $\\beta_{(i,n)}$ correspond to the number of risky securities and the bonds, respectively, for the $i^{\\text{th}}$ path at the $n^{\\text{th}}$ time period. It is to be realized that for a given time period, $1\\leq i \\leq s^n$.\n",
"\n",
"The above-stated representation shows that we can \"repeat\" the variables of the previous level $s$ number of times and then it becomes quite easy to implement the budget constraints in Fusion. One can also extend the above shown method to a higher split index."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"#Each level corresponds to a time step in the event tree. This is a subclass of the 'Tree' class.\n",
"class Level(Tree):\n",
" # l corresponds to the time step; a_old, b_old belong to (l-1) in the tree. \n",
" def __init__(self,l,Tree,a_old,b_old):\n",
" if l==Tree.T :\n",
" #If the current level is the final time period, then all risky securities are considered sold.\n",
" self.a_new = Expr.constTerm([Tree.s**l,Tree.n],0.0)\n",
" #Moreover, the final step's bonds are the wealth W(I_T); later used in the utility function.\n",
" self.b_new = Tree.M.variable('W_T',Tree.s**l)\n",
" else:\n",
" #Risky securities for time step l.\n",
" self.a_new = Tree.M.variable('a_{}'.format(l),[Tree.s**l,Tree.n])\n",
" #Bonds in dollars for time step l.\n",
" self.b_new = Tree.M.variable('b_{}'.format(l),Tree.s**l) \n",
" #Variable for the quadratic cone to implement the absolute value constraint.\n",
" self.t_new = Tree.M.variable('t_{}'.format(l),[Tree.s**l,Tree.n],Domain.greaterThan(0.0))\n",
" Tree.cost = np.repeat(Tree.cost,Tree.s,axis=0)\n",
"\n",
" #Repeating/Stacking the (l-1) level variables to implement the linear budget constraints. \n",
" A = Expr.repeat(a_old,Tree.s,0)\n",
" B = Expr.repeat(b_old,Tree.s,0)\n",
" #The price vector of the previous level, multiplied by the scale factors and then stacked vertically.\n",
" Tree.S_v_t = np.vstack([np.multiply(j,Tree.S_v_t) for j in Tree.S_v])\n",
"\n",
" #Expressions involved in the budget constraints.\n",
" self.bond_sub = Expr.sub(Expr.mul(B,(1+Tree.r)),self.b_new)\n",
" self.secu_sub = Expr.sub(self.a_new,A)\n",
" self.transact = Expr.add(self.secu_sub,Expr.mulElm(Tree.cost,self.t_new))\n",
" self.secu_exp = Expr.mulDiag(self.transact,Tree.S_v_t.transpose())\n",
" \n",
" #The linear budget constraints.\n",
" self.budg_constr = Tree.M.constraint('Budget_{}'.format(l),Expr.sub(self.bond_sub,self.secu_exp),\n",
" Domain.greaterThan(0.0))\n",
" #Quadratic cone for the absolute value requirement in the budget constraints.\n",
" Tree.M.constraint('a{}_abs'.format(l),Expr.stack(2,self.t_new,self.secu_sub),Domain.inQCone())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Now that we have the basic structure of the tree ready, we can easily make the Fusion model as follows..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will represent the paths $I_T$ by their numbers. Therefore, we have (split_index)$^T$ paths. Each Path number represented in the base of the split index will give us a unique id for each path."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"#This function converts the path number from base 10 to base split-index. \n",
"def base_rep(b,i,size):\n",
" n = np.zeros(size)\n",
" for j in range(size):\n",
" n[j] = i%b\n",
" if i//b == 0:\n",
" break\n",
" else:\n",
" i = i//b\n",
" return(np.flip(n))\n",
"\n",
"#This function converts the path number to path id, for all paths.\n",
"def path_id_make(split_index,T):\n",
" path_id = []\n",
" for i in range(split_index**T):\n",
" path_id.append(base_rep(split_index,i,T))\n",
" return(np.asarray(path_id).astype(int))\n",
"\n",
"#This function calculates: path probabilities as well as the price of the underlying stock.\n",
"def path_route_calc(path_id,split_index,q,*args):\n",
" s = np.zeros(split_index)\n",
" for p in path_id:\n",
" s[p] = s[p] + 1\n",
" if args:\n",
" path_S1 = []\n",
" for j in range(split_index):\n",
" path_S1.append((args[0][j])**s[j])\n",
" path_S1T = np.prod(np.asarray(path_S1))\n",
" #Returning price of the stock.\n",
" return(path_S1T)\n",
" else:\n",
" path_prob = np.prod(q**s)\n",
" #Returning probability of the path.\n",
" return(path_prob)\n",
"\n",
"\n",
"#These functions will be used for the calculation of the price-scaling factors (i.e. a,b and c).\n",
"#They follow the discrete approximation of the continuous price vector process.\n",
"def price_vector_z(sigma,eps,dt,q):\n",
" E = np.exp(np.matmul(eps,sigma)*np.sqrt(dt))\n",
" return(np.matmul(q,E))\n",
"\n",
"def price_vector_abc(sigma,eps,dt,q,mu,Z):\n",
" E1 = np.multiply(np.exp(np.matmul(eps,sigma)*np.sqrt(dt)),1/Z)\n",
" E2 = np.repeat([np.exp(mu*dt)],eps.shape[0],axis = 0)\n",
" return(np.multiply(E2,E1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Portfolio without call option:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def Portfolio(port_params,util_type = 'HARA',wrtLog = True):\n",
" T,W_0,S,theta,S_v_0,r,U,q = port_params\n",
" M = Model('PORTFOLIO')\n",
" if wrtLog:\n",
" M.setLogHandler(sys.stdout)\n",
" #Instantiating the Tree class. \n",
" Tree_1 = Tree(M,T,W_0,S,theta,S_v_0,r,U)\n",
" #Making all the levels.\n",
" Tree_1.level_make()\n",
" #Making the conic constraints for the utility function.\n",
" H = Tree_1.util_dispatch[util_type](Tree_1.levels[T-1].b_new)\n",
" #Making the path ids (see function def).\n",
" path_ids = path_id_make(Tree_1.s,T)\n",
" #Calculating the path probabilities (see function def).\n",
" path_probs = [path_route_calc(p,Tree_1.s,q) for p in path_ids]\n",
" #Objective function\n",
" Obj = Expr.dot(H,path_probs)\n",
" M.objective('PORTFOLIO_OBJ',ObjectiveSense.Maximize,Obj)\n",
" M.solve()\n",
" #The maximal expected utility\n",
" utility_W0 = M.primalObjValue()\n",
" ut_time = M.getSolverDoubleInfo('optimizerTime')\n",
" ut_iter = M.getSolverIntInfo('intpntIter')\n",
" return(M,Tree_1,utility_W0,path_ids,Obj,[ut_iter,ut_time])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Portfolio with call option:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This model involves modifying the initial and the final budget constraints, adding a constraint on the new utility, creating a new objective (the option price) and re-optimizing the previous model. Fusion allows re-optimizing (this saves a considerable amount of time, which would otherwise be spent in re-building the model)."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def Option_Portfolio(port_params,K,util_type = 'HARA',writeLog = True,solver_info=False):\n",
" T,W_0,S,theta,S_v_0,r,U,q = port_params\n",
" #We start by solving the problem for the no option case, so that we have U*(W_0).\n",
" M,Tree_1,utility_W0,path_ids,Obj,obj_info = Portfolio(port_params,util_type = util_type,wrtLog = writeLog)\n",
" #This is the price of the underlying stock at the Time horizon (for each path, i.e. s**T paths)\n",
" path_Svs = np.asarray([path_route_calc(p,S.shape[0],q,S[:,0]) for p in path_ids])\n",
" #max((Stock_price(I_T) - K),0)\n",
" call_profit = ((path_Svs - K) + abs(path_Svs - K))/2 \n",
" #Adding a new variable to the model: the reservation purchase price of the option.\n",
" #Note: the call option is bought only on one underlying security.\n",
" Call = M.variable('Call_Price',1,Domain.inRange(0.0,S_v_0[0]))\n",
" #Updating the root constraint, making W_0 into (W_0 - C).\n",
" Tree_1.root_budget.update(Expr.neg(Call),Call)\n",
" #Updating the budget constraints in the final time-period by adding max((Stock_price(I_T) - K),0).\n",
" Tree_1.levels[T-1].budg_constr.update(call_profit.tolist())\n",
" #Now, we add the constraint E[U(W_T)] >= U*(W_0), to our model.\n",
" M.constraint('E_geq_U',Expr.sub(Obj,utility_W0),Domain.greaterThan(0.0))\n",
" #Objective function: given by the reservation purchase price itself\n",
" M.objective('Call_price_OBJECTIVE',ObjectiveSense.Maximize,Call)\n",
" M.solve()\n",
" if solver_info:\n",
" price_time = M.getSolverDoubleInfo('optimizerTime')\n",
" price_iter = M.getSolverIntInfo('intpntIter')\n",
" n_cons = M.getSolverIntInfo('optNumcon')\n",
" n_vars = M.getSolverIntInfo('optNumvar')\n",
" return(M.primalObjValue(),n_cons,n_vars,obj_info[0],obj_info[1],price_iter,price_time)\n",
" else:\n",
" return(M.primalObjValue())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, to demonstrate the model, we solve it for two simple cases. The values taken below are the same as mentioned in the paper (Andersen et. al.). "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Problem\n",
" Name : PORTFOLIO \n",
" Objective sense : max \n",
" Type : CONIC (conic optimization problem)\n",
" Constraints : 2062 \n",
" Cones : 606 \n",
" Scalar variables : 2547 \n",
" Matrix variables : 0 \n",
" Integer variables : 0 \n",
"\n",
"Optimizer started.\n",
"Presolve started.\n",
"Linear dependency checker started.\n",
"Linear dependency checker terminated.\n",
"Eliminator started.\n",
"Freed constraints in eliminator : 152\n",
"Eliminator terminated.\n",
"Eliminator started.\n",
"Freed constraints in eliminator : 0\n",
"Eliminator terminated.\n",
"Eliminator - tries : 2 time : 0.00 \n",
"Lin. dep. - tries : 1 time : 0.00 \n",
"Lin. dep. - number : 0 \n",
"Presolve terminated. Time: 0.01 \n",
"Problem\n",
" Name : PORTFOLIO \n",
" Objective sense : max \n",
" Type : CONIC (conic optimization problem)\n",
" Constraints : 2062 \n",
" Cones : 606 \n",
" Scalar variables : 2547 \n",
" Matrix variables : 0 \n",
" Integer variables : 0 \n",
"\n",
"Optimizer - threads : 20 \n",
"Optimizer - solved problem : the primal \n",
"Optimizer - Constraints : 494\n",
"Optimizer - Cones : 607\n",
"Optimizer - Scalar variables : 1818 conic : 1465 \n",
"Optimizer - Semi-definite variables: 0 scalarized : 0 \n",
"Factor - setup time : 0.00 dense det. time : 0.00 \n",
"Factor - ML order time : 0.00 GP order time : 0.00 \n",
"Factor - nonzeros before factor : 6073 after factor : 6073 \n",
"Factor - dense dim. : 0 flops : 2.67e+05 \n",
"ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME \n",
"0 1.0e+00 1.1e+00 3.2e+02 0.00e+00 0.000000000e+00 3.168332369e+02 1.0e+00 0.04 \n",
"1 3.3e-01 3.8e-01 4.9e+01 1.90e+00 6.579127319e-01 6.906592410e+01 3.3e-01 0.06 \n",
"2 3.0e-02 3.4e-02 1.2e+00 1.55e+00 7.999129416e-01 5.499606733e+00 3.0e-02 0.06 \n",
"3 5.3e-03 6.0e-03 1.1e-01 8.94e-01 1.503700941e+00 2.437452983e+00 5.3e-03 0.07 \n",
"4 9.2e-04 1.1e-03 1.0e-02 6.83e-01 2.140975891e+00 2.341291055e+00 9.2e-04 0.07 \n",
"5 1.1e-04 1.2e-04 5.2e-04 7.13e-01 2.517001446e+00 2.544764609e+00 1.1e-04 0.08 \n",
"6 2.7e-05 3.0e-05 6.3e-05 9.50e-01 2.568634524e+00 2.575638817e+00 2.7e-05 0.08 \n",
"7 2.0e-05 2.3e-05 4.2e-05 8.04e-01 2.573149607e+00 2.578786129e+00 2.0e-05 0.09 \n",
"8 1.1e-05 1.3e-05 1.8e-05 7.55e-01 2.580098187e+00 2.583562469e+00 1.1e-05 0.10 \n",
"9 3.6e-06 4.1e-06 3.8e-06 7.46e-01 2.587301575e+00 2.588645267e+00 3.6e-06 0.10 \n",
"10 1.8e-06 2.0e-06 1.4e-06 5.83e-01 2.590268239e+00 2.591037107e+00 1.8e-06 0.11 \n",
"11 7.4e-07 8.5e-07 3.7e-07 7.71e-01 2.592230984e+00 2.592572571e+00 7.4e-07 0.11 \n",
"12 1.0e-07 1.2e-07 1.8e-08 9.09e-01 2.593476774e+00 2.593525346e+00 1.0e-07 0.11 \n",
"13 8.3e-10 9.5e-10 1.2e-11 9.94e-01 2.593674050e+00 2.593674447e+00 8.3e-10 0.12 \n",
"14 8.3e-09 1.1e-10 5.1e-13 1.00e+00 2.593675395e+00 2.593675443e+00 1.0e-10 0.12 \n",
"15 7.8e-09 1.1e-10 4.7e-13 9.83e-01 2.593675406e+00 2.593675452e+00 9.5e-11 0.12 \n",
"16 8.4e-09 9.4e-11 3.8e-13 1.00e+00 2.593675430e+00 2.593675469e+00 8.3e-11 0.13 \n",
"17 7.6e-09 7.3e-11 2.6e-13 1.00e+00 2.593675464e+00 2.593675495e+00 6.4e-11 0.15 \n",
"Optimizer terminated. Time: 0.15 \n",
"\n",
"\n",
"Interior-point solution summary\n",
" Problem status : PRIMAL_AND_DUAL_FEASIBLE\n",
" Solution status : OPTIMAL\n",
" Primal. obj: 2.5936754639e+00 nrm: 6e+00 Viol. con: 4e-08 var: 0e+00 cones: 0e+00 \n",
" Dual. obj: 2.5936754946e+00 nrm: 3e+00 Viol. con: 2e-11 var: 1e-10 cones: 0e+00 \n",
"Optimizer started.\n",
"Optimizer terminated. Time: 0.09 \n",
"\n",
"\n",
"Interior-point solution summary\n",
" Problem status : PRIMAL_AND_DUAL_FEASIBLE\n",
" Solution status : OPTIMAL\n",
" Primal. obj: 5.3868459819e-02 nrm: 6e+00 Viol. con: 3e-06 var: 0e+00 cones: 2e-07 \n",
" Dual. obj: 5.3868462409e-02 nrm: 4e+00 Viol. con: 1e-12 var: 9e-12 cones: 0e+00 \n",
"\n",
"\n",
"Call option price = 0.0538684598192706\n"
]
}
],
"source": [
"T = 5 #Number of periods\n",
"sigma = np.asarray([0.2]) #Sigma matrix\n",
"q = np.ones(3)/3 #Probabilities\n",
"eps = np.asarray([[np.sqrt(2)],[-1/np.sqrt(2)],[-1/np.sqrt(2)]]) #Eps matrix (see text)\n",
"T_horizon = 1/4 #Time horizon (in years)\n",
"dT = T_horizon/T \n",
"theta = [0.005] #Transaction costs\n",
"S_v_0 = [1.0] #Initial price vector\n",
"r = (1.06**dT) - 1 #Constant interest rate\n",
"W_0 = 1.0 #Initial wealth\n",
"gamma = 0.3 #Gamma parameter for HARA utility\n",
"b = 1 #b parameter for HARA utility (see text)\n",
"c = 0.2 #Absolute risk aversion ARA(W)\n",
"a = b/((1/c) + (W_0/(gamma-1)))\n",
"K = 1 #Strike Price.\n",
"mu = 0.15 #Drift of the stock price\n",
"\n",
"Z = price_vector_z(sigma,eps,dT,q)\n",
"S_coeffs = np.expand_dims(price_vector_abc(sigma,eps,dT,q,mu,Z),axis=1)\n",
" \n",
"util_paras = np.asarray([a,b,gamma])\n",
"input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,util_paras,q]\n",
"\n",
"#Using HARA utility\n",
"print('\\n\\nCall option price = {}'.format(Option_Portfolio(input_pars,K)))"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Problem\n",
" Name : PORTFOLIO \n",
" Objective sense : max \n",
" Type : CONIC (conic optimization problem)\n",
" Constraints : 1819 \n",
" Cones : 606 \n",
" Scalar variables : 2547 \n",
" Matrix variables : 0 \n",
" Integer variables : 0 \n",
"\n",
"Optimizer started.\n",
"Presolve started.\n",
"Linear dependency checker started.\n",
"Linear dependency checker terminated.\n",
"Eliminator started.\n",
"Freed constraints in eliminator : 152\n",
"Eliminator terminated.\n",
"Eliminator started.\n",
"Freed constraints in eliminator : 0\n",
"Eliminator terminated.\n",
"Eliminator - tries : 2 time : 0.00 \n",
"Lin. dep. - tries : 1 time : 0.00 \n",
"Lin. dep. - number : 0 \n",
"Presolve terminated. Time: 0.01 \n",
"Problem\n",
" Name : PORTFOLIO \n",
" Objective sense : max \n",
" Type : CONIC (conic optimization problem)\n",
" Constraints : 1819 \n",
" Cones : 606 \n",
" Scalar variables : 2547 \n",
" Matrix variables : 0 \n",
" Integer variables : 0 \n",
"\n",
"Optimizer - threads : 20 \n",
"Optimizer - solved problem : the primal \n",
"Optimizer - Constraints : 494\n",
"Optimizer - Cones : 607\n",
"Optimizer - Scalar variables : 1818 conic : 1465 \n",
"Optimizer - Semi-definite variables: 0 scalarized : 0 \n",
"Factor - setup time : 0.00 dense det. time : 0.00 \n",
"Factor - ML order time : 0.00 GP order time : 0.00 \n",
"Factor - nonzeros before factor : 6073 after factor : 6073 \n",
"Factor - dense dim. : 0 flops : 2.67e+05 \n",
"ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME \n",
"0 1.0e+00 1.3e+00 2.0e+02 0.00e+00 -1.454638549e+00 2.006397864e+02 1.0e+00 0.01 \n",
"1 5.7e-01 7.2e-01 9.0e+01 6.27e-01 -2.672174116e+00 1.209604619e+02 5.7e-01 0.02 \n",
"2 1.6e-01 2.0e-01 1.1e+01 8.90e-01 -4.662276868e+00 2.902021528e+01 1.6e-01 0.03 \n",
"3 2.5e-02 3.2e-02 3.5e-01 1.48e+00 -4.266144780e+00 -4.202866560e-01 2.5e-02 0.03 \n",
"4 5.9e-03 7.5e-03 3.9e-02 1.78e+00 2.263613352e-01 8.439751686e-01 5.9e-03 0.03 \n",
"5 8.1e-04 1.0e-03 1.9e-03 1.15e+00 8.112108317e-01 8.894400473e-01 8.1e-04 0.03 \n",
"6 2.5e-04 3.1e-04 3.4e-04 1.03e+00 8.873422767e-01 9.111334760e-01 2.5e-04 0.04 \n",
"7 1.4e-04 1.8e-04 1.7e-04 8.14e-01 9.026036581e-01 9.179474679e-01 1.4e-04 0.04 \n",
"8 7.2e-05 9.1e-05 7.5e-05 6.46e-01 9.156954394e-01 9.253195920e-01 7.2e-05 0.04 \n",
"9 3.1e-05 4.0e-05 2.6e-05 5.95e-01 9.273205811e-01 9.324555987e-01 3.1e-05 0.04 \n",
"10 1.4e-05 1.7e-05 1.0e-05 4.18e-01 9.368898705e-01 9.399326786e-01 1.4e-05 0.04 \n",
"11 9.5e-06 1.2e-05 6.3e-06 5.92e-01 9.404745780e-01 9.427436466e-01 9.5e-06 0.05 \n",
"12 6.1e-06 7.8e-06 3.4e-06 6.52e-01 9.438858106e-01 9.454515285e-01 6.1e-06 0.05 \n",
"13 2.5e-06 3.2e-06 9.8e-07 7.56e-01 9.480435837e-01 9.487414991e-01 2.5e-06 0.05 \n",
"14 4.8e-07 6.1e-07 8.3e-08 8.98e-01 9.508842290e-01 9.510215048e-01 4.8e-07 0.05 \n",
"15 3.8e-08 4.8e-08 1.9e-09 9.86e-01 9.515180580e-01 9.515289900e-01 3.8e-08 0.05 \n",
"16 1.9e-09 2.4e-09 2.2e-11 9.99e-01 9.515724494e-01 9.515729966e-01 1.9e-09 0.05 \n",
"17 4.1e-09 1.8e-10 4.4e-13 1.00e+00 9.515750903e-01 9.515751312e-01 1.4e-10 0.06 \n",
"Optimizer terminated. Time: 0.06 \n",
"\n",
"\n",
"Interior-point solution summary\n",
" Problem status : PRIMAL_AND_DUAL_FEASIBLE\n",
" Solution status : OPTIMAL\n",
" Primal. obj: 9.5157509030e-01 nrm: 6e+00 Viol. con: 8e-09 var: 0e+00 cones: 0e+00 \n",
" Dual. obj: 9.5157513120e-01 nrm: 1e+00 Viol. con: 1e-16 var: 3e-10 cones: 0e+00 \n",
"Optimizer started.\n",
"Optimizer terminated. Time: 0.04 \n",
"\n",
"\n",
"Interior-point solution summary\n",
" Problem status : PRIMAL_AND_DUAL_FEASIBLE\n",
" Solution status : OPTIMAL\n",
" Primal. obj: 5.3526766778e-02 nrm: 6e+00 Viol. con: 7e-08 var: 0e+00 cones: 9e-09 \n",
" Dual. obj: 5.3526769271e-02 nrm: 1e+00 Viol. con: 7e-17 var: 2e-11 cones: 0e+00 \n",
"\n",
"\n",
"Call option price = 0.053526766777519254\n"
]
}
],
"source": [
"#Using Exponential utility (everything else is the same as above)\n",
"input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]\n",
"print('\\n\\nCall option price = {}'.format(Option_Portfolio(input_pars,K,util_type='EXP')))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Test-cases:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Reservation price sensitivity for the choice of utility function:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(a.) Dependency of the reservation purchase price on $\\gamma$, $ARA(W_0)$ and $\\overline{T}$. (Consider Table 3, Andersen et. al.)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ARA = 0.2 ; Time Horizon = 1/4 year\n",
"gamma C \n",
"-4.0 0.0535872\n",
"-2.0 0.0536078\n",
"-0.9 0.0536516\n",
"-0.3 0.0537167\n",
" 0.3 0.0538685\n",
" 0.6 0.0541304\n",
" exp 0.0535268\n"
]
}
],
"source": [
"gamma = [-4.0,-2.0,-0.9,-0.3,0.3,0.6] #Gamma parameter for HARA utility\n",
"\n",
"print('ARA = 0.2 ; Time Horizon = 1/4 year')\n",
"print('{0:^6} {1:^9}'.format('gamma','C'))\n",
"for g in gamma:\n",
" a = b/((1/c) + (W_0/(g-1)))\n",
" util_paras = np.asarray([a,b,g])\n",
" input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,util_paras,q]\n",
" #Using HARA utility\n",
" print('{0:^ 5.1f} {1:^9.7f}'.format(g,Option_Portfolio(input_pars,K,writeLog=False)))\n",
"input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]\n",
"#Using Exponential utility\n",
"print('{0:^5} {1:^9.7f}'.format('exp',Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ARA = 1.0 ; Time Horizon = 1/4 year\n",
"gamma C \n",
"-4.0 0.0533297\n",
"-2.0 0.0533626\n",
"-0.9 0.0534029\n",
"-0.3 0.0534659\n",
" exp 0.0532800\n"
]
}
],
"source": [
"c = 1.0 #Absolute risk aversion ARA(W)\n",
"\n",
"print('ARA = 1.0 ; Time Horizon = 1/4 year')\n",
"print('{0:^6} {1:^9}'.format('gamma','C'))\n",
"for g in gamma[0:4]:\n",
" a = b/((1/c) + (W_0/(g-1)))\n",
" util_paras = np.asarray([a,b,g])\n",
" input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,util_paras,q]\n",
" #Using HARA utility\n",
" print('{0:^ 5.1f} {1:^9.7f}'.format(g,Option_Portfolio(input_pars,K,writeLog=False)))\n",
"input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]\n",
"#Using Exponential utility\n",
"print('{0:^5} {1:^9.7f}'.format('exp',Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ARA = 0.2 ; Time Horizon = 9 years\n",
"gamma C \n",
"-4.0 0.4182626\n",
"-2.0 0.4181730\n",
"-0.9 0.4181584\n",
"-0.3 0.4182109\n",
" 0.3 0.4181520\n",
" 0.6 0.4128452\n",
" exp 0.4181519\n"
]
}
],
"source": [
"T_horizon = 9 #Time horizon (in years)\n",
"dT = T_horizon/T \n",
"Z = price_vector_z(sigma,eps,dT,q)\n",
"S_coeffs = np.expand_dims(price_vector_abc(sigma,eps,dT,q,mu,Z),axis=1) \n",
"r = (1.06**dT) - 1 #Constant interest rate\n",
"c = 0.2 #Absolute risk aversion ARA(W)\n",
"\n",
"print('ARA = 0.2 ; Time Horizon = 9 years')\n",
"print('{0:^6} {1:^9}'.format('gamma','C'))\n",
"for g in gamma:\n",
" a = b/((1/c) + (W_0/(g-1)))\n",
" util_paras = np.asarray([a,b,g])\n",
" input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,util_paras,q]\n",
" #Using HARA utility\n",
" print('{0:^ 5.1f} {1:^9.7f}'.format(g,Option_Portfolio(input_pars,K,writeLog=False)))\n",
"input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]\n",
"#Using Exponential utility\n",
"print('{0:^5} {1:^9.7f}'.format('exp',Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ARA = 1.0 ; Time Horizon = 9 years\n",
"gamma C \n",
"-4.0 0.4181606\n",
"-2.0 0.4181575\n",
"-0.9 0.4181590\n",
"-0.3 0.4181670\n",
" exp 0.4176740\n"
]
}
],
"source": [
"c = 1.0 #Absolute risk aversion ARA(W)\n",
"\n",
"print('ARA = 1.0 ; Time Horizon = 9 years')\n",
"print('{0:^6} {1:^9}'.format('gamma','C'))\n",
"for g in gamma[0:4]:\n",
" a = b/((1/c) + (W_0/(g-1)))\n",
" util_paras = np.asarray([a,b,g])\n",
" input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,util_paras,q]\n",
" #Using HARA utility\n",
" print('{0:^ 5.1f} {1:^9.7f}'.format(g,Option_Portfolio(input_pars,K,writeLog=False)))\n",
"input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]\n",
"#Using Exponential utility\n",
"print('{0:^5} {1:^9.7f}'.format('exp',Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(b.) Dependency of the reservation purchase price on the initial level of Relative Risk Aversion, $RRA(W_0)$ (Consider Table 4, Andersen et. al.)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ARA = 0.2 ; Time Horizon = 1/4 year\n",
" W_0 RRA C \n",
" 1 0.2 0.0536405\n",
" 4 0.8 0.0536451\n",
" 8 1.6 0.0540775\n",
" exp 0.0535268\n"
]
}
],
"source": [
"T = 5 #Number of periods\n",
"T_horizon = 1/4 #Time horizon (in years)\n",
"dT = T_horizon/T \n",
"r = (1.06**dT) - 1 #Constant interest rate\n",
"W_0 = [1.0,4.0,8.0] #Initial wealth\n",
"gamma = -1.0 #Gamma parameter for HARA utility\n",
"b = 1 #b parameter for HARA utility (see text)\n",
"c = 0.2 #Absolute risk aversion ARA(W)\n",
"Z = price_vector_z(sigma,eps,dT,q)\n",
"S_coeffs = np.expand_dims(price_vector_abc(sigma,eps,dT,q,mu,Z),axis=1)\n",
"\n",
"print('ARA = 0.2 ; Time Horizon = 1/4 year')\n",
"print('{0:^5} {1:^5} {2:^9}'.format('W_0','RRA','C'))\n",
"for w0 in W_0:\n",
" a = b/((1/c) + (w0/(gamma-1)))\n",
" util_paras = np.asarray([a,b,gamma])\n",
" input_pars = [T,w0,S_coeffs,theta,S_v_0,r,util_paras,q]\n",
" #Using HARA utility\n",
" print('{0:^5.0f} {1:^5.1f} {2:^9.7f}'.format(w0,c*w0,Option_Portfolio(input_pars,K,writeLog=False)))\n",
"input_pars = [T,1.0,S_coeffs,theta,S_v_0,r,c,q]\n",
"#Using Exponential utility\n",
"print('{0:^5} {1:^5} {2:^9.7f}'.format('exp','',Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ARA = 0.2 ; Time Horizon = 9 years\n",
" W_0 RRA C \n",
" 1 0.2 0.4181970\n",
" 4 0.8 0.4182284\n",
" 8 1.6 0.4183691\n",
" exp 0.4181519\n"
]
}
],
"source": [
"T_horizon = 9 #Time horizon (in years)\n",
"dT = T_horizon/T\n",
"r = (1.06**dT) - 1 #Constant interest rate\n",
"Z = price_vector_z(sigma,eps,dT,q)\n",
"S_coeffs = np.expand_dims(price_vector_abc(sigma,eps,dT,q,mu,Z),axis=1)\n",
"\n",
"print('ARA = 0.2 ; Time Horizon = 9 years')\n",
"print('{0:^5} {1:^5} {2:^9}'.format('W_0','RRA','C'))\n",
"for w0 in W_0:\n",
" a = b/((1/c) + (w0/(gamma-1)))\n",
" util_paras = np.asarray([a,b,gamma])\n",
" input_pars = [T,w0,S_coeffs,theta,S_v_0,r,util_paras,q]\n",
" #Using HARA utility\n",
" print('{0:^5.0f} {1:^5.1f} {2:^9.7f}'.format(w0,c*w0,Option_Portfolio(input_pars,K,writeLog=False)))\n",
"input_pars = [T,1.0,S_coeffs,theta,S_v_0,r,c,q]\n",
"\n",
"print('{0:^5} {1:^5} {2:^9.7f}'.format('exp','',Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(c.) Convergence of the reservation purchase price. (Consider Table 5, Andersen et. al.)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" -exp(-0.7W) (W^0.3)/0.3 \n",
" T C_exp C_pow CRR \n",
" 1 0.0513874752 0.0513875149 0.0429986634\n",
" 2 0.0566847213 0.0568627676 0.0498155386\n",
" 3 0.0533733615 0.0536936823 0.0471563891\n",
" 4 0.0491337501 0.0494956789 0.0429546099\n",
" 5 0.0533720359 0.0537069769 0.0471967162\n",
" 6 0.0533763586 0.0537373351 0.0473670422\n",
" 7 0.0509534861 0.0513439146 0.0450257209\n",
" 8 0.0524868169 0.0528617415 0.0465173011\n",
" 9 0.0531665811 0.0535465362 0.0472606223\n",
"10 0.0519294335 0.0523270335 0.0460868513\n",
"11 0.0520640889 0.0524579600 0.0461999995\n"
]
}
],
"source": [
"T_horizon = 1/4 #Time horizon (in years)\n",
"W_0 = 1.0 #Initial wealth\n",
"gamma = 0.3 #Gamma parameter for HARA utility\n",
"b = 0.0 #b parameter for HARA utility (see text)\n",
"c = 0.7 #Absolute risk aversion ARA(W)\n",
"a = 1 - gamma\n",
"util_paras = np.asarray([a,b,gamma])\n",
"\n",
"C_trio = []\n",
"print('{0:^3} {1:12} {2:12} {3:^12}'.format(' ','-exp(-0.7W)','(W^0.3)/0.3',' '))\n",
"print('{0:^3} {1:^12} {2:^12} {3:^12}'.format('T','C_exp','C_pow','CRR'))\n",
"for T in np.arange(1,12):\n",
" dT = T_horizon/T\n",
" r = (1.06**dT) - 1\n",
" Z = price_vector_z(sigma,eps,dT,q)\n",
" S_coeffs = np.expand_dims(price_vector_abc(sigma,eps,dT,q,mu,Z),axis=1)\n",
" #C_pow : Power utility function.\n",
" input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,util_paras,q]\n",
" C_pow = Option_Portfolio(input_pars,K,writeLog=False)\n",
" #CRR : Call price in a frictionless market \n",
" input_pars = [T,W_0,S_coeffs,[0.0],S_v_0,r,util_paras,q]\n",
" CRR = Option_Portfolio(input_pars,K,writeLog=False)\n",
" #C_exp : Exponential utility function.\n",
" input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]\n",
" C_exp = Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)\n",
" print('{0:^3d} {1:^12.10f} {2:^12.10f} {3:^12.10f}'.format(T,C_exp,C_pow,CRR))\n",
" C_trio.append([C_exp,C_pow,CRR])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(d.) Reservation purchase price dependence on absolute risk aversion (Exponential utility). (Consider Table 6, Andersen et. al.)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ARA(W_0) C[ARA(W_0)] \n",
" 0.10 0.0533767332\n",
" 0.20 0.0533414293\n",
" 0.40 0.0532712772\n",
" 0.80 0.0531318449\n",
" 1.00 0.0530625844\n",
" 2.00 0.0527035845\n",
" 4.00 0.0503777401\n"
]
}
],
"source": [
"T = 9 #Number of periods\n",
"dT = T_horizon/T \n",
"theta = [0.005] #Transaction costs\n",
"r = (1.06**dT) - 1 #Constant interest rate\n",
"c = [0.1,0.2,0.4,0.8,1.0,2.0,4.0] #Absolute risk aversion ARA(W)\n",
"Z = price_vector_z(sigma,eps,dT,q)\n",
"S_coeffs = np.expand_dims(price_vector_abc(sigma,eps,dT,q,mu,Z),axis=1) \n",
"\n",
"print('{0:^8} {1:^12}'.format('ARA(W_0)','C[ARA(W_0)]'))\n",
"C_ara = []\n",
"for kappa in c:\n",
" input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,kappa,q]\n",
" C_exp = Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)\n",
" print('{0:^8.2f} {1:^12.10f}'.format(kappa,C_exp))\n",
" C_ara.append(C_exp)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(e.) Reservation purchase price dependence on absolute risk aversion (Power utility). (Consider Table 7, Andersen et. al.)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" eta C[ARA(W_0)] \n",
"0.10 0.0550906447\n",
"0.20 0.0545760010\n",
"0.30 0.0542082969\n",
"0.40 0.0539549889\n",
"0.50 0.0537753439\n",
"0.60 0.0536467241\n",
"0.70 0.0535465036\n",
"0.80 0.0534635107\n",
"0.90 0.0533917250\n"
]
}
],
"source": [
"eta = np.linspace(0.1,0.9,9) #Risk Aversion parameter\n",
"a = 1\n",
"b = 0.0\n",
"\n",
"print('{0:^5} {1:^12}'.format('eta','C[ARA(W_0)]'))\n",
"C_eta = []\n",
"for n in eta:\n",
" gamma = 1-n\n",
" util_paras = np.asarray([a,b,gamma])\n",
" input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,util_paras,q]\n",
" C_pow = Option_Portfolio(input_pars,K,writeLog=False)\n",
" print('{0:^5.2f} {1:^12.10f}'.format(n,C_pow))\n",
" C_eta.append(C_pow)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(f.) Reservation purchase price dependence on initial wealth (Power utility with $\\eta = 0.7$). (Consider Table 8, Andersen et. al.)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" W_0 C_pow(W_0) \n",
"0.25 0.0525390204\n",
"0.50 0.0532969822\n",
"1.00 0.0535464882\n",
"2.00 0.0536744163\n",
"4.00 0.0537396292\n",
"8.00 0.0537739025\n"
]
}
],
"source": [
"eta = 0.7 #Risk Aversion parameter\n",
"W_0 = [0.25,0.50,1.0,2.0,4.0,8.0] #Initial wealth\n",
"gamma = 1-eta\n",
"util_paras = np.asarray([a,b,gamma])\n",
"\n",
"print('{0:^5} {1:^12}'.format('W_0','C_pow(W_0)'))\n",
"C_w0 = []\n",
"for w0 in W_0:\n",
" input_pars = [T,w0,S_coeffs,theta,S_v_0,r,util_paras,q]\n",
" C_pow = Option_Portfolio(input_pars,K,writeLog=False)\n",
" print('{0:^5.2f} {1:^12.10f}'.format(w0,C_pow))\n",
" C_w0.append(C_pow)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The effect of diversification opportunities on the reservation purchase price."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The model presented above is fairly general, and therefore we can easily extend its application to the case where trade takes place in multiple securities. For illustration, we present the case where there are two risky securities, following Andersen et. al.. For this example, we have:\n",
"\n",
"$$\n",
"\\sigma^{T} = \\begin{pmatrix}\n",
"\\sigma_{11}\\, ,\\,\\sigma_{21} \\\\\n",
"\\sigma_{12}\\, ,\\,\\sigma_{22} \\\\\n",
"\\end{pmatrix} = \\begin{pmatrix}\n",
"\\sigma_{11}\\, ,\\, 0.2\\rho\\\\\n",
"0.0\\, , \\, \\sqrt{(0.2)^2 - \\sigma_{21}^2} \\\\\n",
"\\end{pmatrix}\n",
"$$\n",
"\n",
"where $\\rho$ is the correlation between logarithms of the two securities. The $\\epsilon$ matrix is given by, \n",
"\n",
"$$\\epsilon = \\begin{pmatrix}\n",
"\\epsilon_1 (\\omega_1),\\epsilon_2 (\\omega_1) \\\\\n",
"\\epsilon_1 (\\omega_2),\\epsilon_2 (\\omega_2) \\\\\n",
"\\epsilon_1 (\\omega_3),\\epsilon_2 (\\omega_3)\n",
"\\end{pmatrix} = \\begin{pmatrix}\n",
"\\sqrt{2}\\,\\,,\\,\\,0 \\\\\n",
"-1/\\sqrt{2},\\sqrt{3/2} \\\\\n",
"-1/\\sqrt{2},\\sqrt{3/2}\n",
"\\end{pmatrix}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Moreover, we will consider different values for the transaction costs. There are three situations considered below: completely frictionless market, friction in one security, friction in both securities.\n",
"\n",
"Note: In the code, we use the transpose of the sigma matrix, hence the shapes of the arrays that represent $\\sigma$ and $\\epsilon$ must be maintained similar to what is shown in the code below."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(g.) Reservation purchase price of a call option for T=9 as a function of $\\rho$ and $\\theta_i$ in the presence of two risky securities. (Consider Table 10, Andersen et. al.)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"theta_1 = 0.000 ; theta_2 = 0.000\n",
" rho C(S1,S2) \n",
" 1.0 0.0472606153\n",
" 0.5 0.0472606445\n",
" 0.0 0.0472605948\n",
"-0.5 0.0472605957\n",
"-0.9 0.0472606051\n"
]
}
],
"source": [
"#No transaction costs (Friction-less market)\n",
"W_0 = 1.0 #Initial wealth\n",
"gamma = 0.3 #Gamma parameter for HARA utility\n",
"b = 1 #b parameter for HARA utility (see text)\n",
"c = 0.2 #Absolute risk aversion ARA(W)\n",
"a = b/((1/c) + (W_0/(gamma-1)))\n",
"util_paras = np.asarray([a,b,gamma])\n",
"K = 1 #Strike Price.\n",
"\n",
"#Parameters for the 2 risky securities under consideration:\n",
"#The sigma matrix entries for the first security are same as before...\n",
"sig11 = 0.2 #sigma_{1,1}\n",
"sig12 = 0.0 #sigma_{1,2}\n",
"#Epsilon matrix (Note the extra column for the second risky security)\n",
"eps = np.asarray([[np.sqrt(2),0],[-1/np.sqrt(2),np.sqrt(3/2)],[-1/np.sqrt(2),-np.sqrt(3/2)]])\n",
"S_v_0 = [1.0,1.0] #Initial price vector\n",
"mu = np.asarray([0.15,0.15]) #Drift of the stock price\n",
"rho = [1,0.5,0.0,-0.5,-0.9] #Correlation\n",
"\n",
"theta = [0.0,0.0] #Transaction costs\n",
"print('theta_1 = {0:5.3f} ; theta_2 = {1:5.3f}'.format(theta[0],theta[1]))\n",
"print('{0:^5} {1:^12}'.format('rho','C(S1,S2)'))\n",
"for p in rho:\n",
" sig21 = p*0.2\n",
" #Sigma matrix\n",
" sigma = np.asarray([[sig11,sig21],[sig12,np.sqrt((0.2**2) - (sig21**2))]])\n",
" Z = price_vector_z(sigma,eps,dT,q)\n",
" S_coeffs = price_vector_abc(sigma,eps,dT,q,mu,Z)\n",
" input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]\n",
" print('{0:^ 4.1f} {1:^12.10f}'.format(p,Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"theta_1 = 0.005 ; theta_2 = 0.000\n",
" rho C(S1,S2) \n",
" 1.00 0.0472606924\n",
" 0.50 0.0535486788\n",
" 0.00 0.0533472924\n",
"-0.50 0.0531698739\n",
"-0.90 0.0530348364\n"
]
}
],
"source": [
"theta = [0.005,0.0] #Transaction costs\n",
"print('theta_1 = {0:5.3f} ; theta_2 = {1:5.3f}'.format(theta[0],theta[1]))\n",
"print('{0:^5} {1:^12}'.format('rho','C(S1,S2)'))\n",
"for p in rho:\n",
" sig21 = p*0.2\n",
" #Sigma matrix\n",
" sigma = np.asarray([[sig11,sig21],[sig12,np.sqrt((0.2**2) - (sig21**2))]])\n",
" Z = price_vector_z(sigma,eps,dT,q)\n",
" S_coeffs = price_vector_abc(sigma,eps,dT,q,mu,Z)\n",
" input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]\n",
" print('{0:^ 4.2f} {1:^12.10f}'.format(p,Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"theta_1 = 0.005 ; theta_2 = 0.005\n",
" rho C(S1,S2) \n",
" 1.00 0.0533413151\n",
" 0.50 0.0533755301\n",
" 0.00 0.0533420533\n",
"-0.50 0.0530762258\n",
"-0.90 0.0520002256\n"
]
}
],
"source": [
"theta = [0.005,0.005] #Transaction costs\n",
"print('theta_1 = {0:5.3f} ; theta_2 = {1:5.3f}'.format(theta[0],theta[1]))\n",
"print('{0:^5} {1:^12}'.format('rho','C(S1,S2)'))\n",
"for p in rho:\n",
" sig21 = p*0.2\n",
" #Sigma matrix\n",
" sigma = np.asarray([[sig11,sig21],[sig12,np.sqrt((0.2**2) - (sig21**2))]])\n",
" Z = price_vector_z(sigma,eps,dT,q)\n",
" S_coeffs = price_vector_abc(sigma,eps,dT,q,mu,Z)\n",
" input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]\n",
" print('{0:^ 4.2f} {1:^12.10f}'.format(p,Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(h.) Reservation purchase price of a call option for T=9, as a function of $\\theta_i$, in the case of two risky securities; $\\rho = -0.9$. (Consider Table 11, Andersen et. al.)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"theta C(S1,S2) \n",
"0.0008 0.0477326668\n",
"0.0016 0.0483956225\n",
"0.0030 0.0497638572\n",
"0.0060 0.0532547315\n",
"0.0100 0.0595919990\n",
"0.0200 0.0618326015\n",
"0.0500 0.0618326729\n",
"0.1000 0.0618325714\n"
]
}
],
"source": [
"rho = -0.9 #Correlation\n",
"sig21 = rho*0.2\n",
"sigma = np.asarray([[sig11,sig21],[sig12,np.sqrt((0.2**2) - (sig21**2))]])\n",
"Z = price_vector_z(sigma,eps,dT,q)\n",
"S_coeffs = price_vector_abc(sigma,eps,dT,q,mu,Z)\n",
"\n",
"theta_list = [0.0008,0.0016,0.003,0.006,0.01,0.02,0.05,0.1]\n",
"print('{0:^5} {1:^12}'.format('theta','C(S1,S2)'))\n",
"for theta_i in theta_list:\n",
" #We take the transaction costs for both securities to be equal.\n",
" theta = [theta_i,theta_i]\n",
" input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]\n",
" print('{0:^5.4f} {1:^12.10f}'.format(theta_i,Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(i.) Computational efficiency for different values of T (time periods). (Consider Table 12, Andersen et. al.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The columns of the output table in the following cell denote: T (time periods), cons. (number of constraints), vars. (number of variables), it. (number of iterations for the no-option model as well as the option-price model), time. (time taken, in seconds for both models)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" No-option Option-price\n",
" T cons. vars. it. time it. time \n",
" 1 26 38 8 0.002 9 0.002\n",
" 2 89 131 7 0.003 8 0.003\n",
" 3 278 410 8 0.005 12 0.006\n",
" 4 845 1247 17 0.023 19 0.022\n",
" 5 2546 3758 20 0.062 21 0.056\n",
" 6 7649 11291 24 0.169 25 0.164\n",
" 7 22958 33890 26 0.672 25 0.499\n",
" 8 68885 101687 30 2.594 28 1.648\n",
" 9 206666 305078 46 22.571 34 8.420\n",
"10 620009 915251 47 50.372 40 29.011\n"
]
}
],
"source": [
"theta = [0.005,0.005]\n",
"print('{0:^2} {1:^8} {2:^10} {3:^11} {4:^11}'.format('','','','No-option','Option-price'))\n",
"print('{0:>2} {1:>8} {2:>10} {3:^3} {4:^6} {5:^3} {6:^6}'.\n",
" format('T','cons.','vars.','it.','time','it.','time'))\n",
"total_info_arr = []\n",
"for T in range(1,10):\n",
" try:\n",
" T_horizon = 1/4\n",
" dt = T_horizon/T\n",
" r = (1.06**dT) - 1\n",
" Z = price_vector_z(sigma,eps,dT,q)\n",
" S_coeffs = price_vector_abc(sigma,eps,dT,q,mu,Z)\n",
" input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]\n",
" call,n_cons,n_var,ut_it,ut_time,price_it,price_time = Option_Portfolio(input_pars,K,util_type='EXP',\n",
" writeLog=False,solver_info=True)\n",
" print('{0:>2d} {1:>8d} {2:>10d} {3:>3d} {4:>6.3f} {5:>3d} {6:>6.3f}'.\n",
" format(T,n_cons,n_var,ut_it,ut_time,price_it,price_time))\n",
" total_info_arr.append([call,n_cons,n_var,ut_it,ut_time,price_it,price_time])\n",
" except MemoryError as e:\n",
" print(e)\n",
" break\n",
" except SolutionError as s:\n",
" print(s)\n",
" except OptimizerError as e:\n",
" print(e)\n",
" break"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The calculations shown above were performed on a desktop with 15.6 GB of RAM and an Intel$^\\circledR$ Core$^{\\text{TM}}$ i7-6770HQ CPU @ 2.6 GHz $\\times$ 8."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
This work is licensed under a Creative Commons Attribution 4.0 International License. The **MOSEK** logo and name are trademarks of Mosek ApS. The code is provided as-is. Compatibility with future release of **MOSEK** or the `Fusion API` are not guaranteed. For more information contact our [support](mailto:support@mosek.com). "
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}