#!/usr/bin/env python
# coding: utf-8

# ![MOSEK ApS](https://www.mosek.com/static/images/branding/webgraphmoseklogocolor.png )

# # Utility based option pricing with transaction costs and diversification

# In this notebook, we present the <i>Fusion</i> 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. 

# ## The Model

# 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: 
# 
# $$(S_{1,n},S_{2,n},..,S_{m,n}) = \begin{cases}
# (a_1 S_{1,n-1}, a_2 S_{2,n-1},..,a_m S_{m,n-1}), \,\, \text{with probability } q_1, \\
# (b_1 S_{1,n-1}, b_2 S_{2,n-1},..,b_m S_{m,n-1}), \,\, \text{with probability } q_2, \\
# (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)
# \end{cases}
# $$
# 
# 
# 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.

# 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 <i>'budget constraint'</i> on the bonds an investor can buy in the next period is given by:
# 
# 
# $$\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)|]$$
# 
# 
# 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:
# 
# $$\beta_0 \leq W_0 - \sum_{i=1}^{m} S_{i,0}[\alpha_{i,0} + |\alpha_{i,0}|]$$
# 
# 
# 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:
# 
# $$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)|]$$

# ### 1.) Portfolio choice

# 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:
# 
# $$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))$$
# 
# 
# 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:
# 
# 
# $$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)]$$
# 
# 
# $$\text{subject to:  } \beta_0 \leq W_0 - \sum_{i=1}^{m} S_{i,0}[\alpha_{i,0} + |\alpha_{i,0}|],$$
# $$\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)|],$$
# $$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$$

# ### 2.) Price vector process

# The continuous price process that we intended to approximate is:
# 
# $$\frac{d\tilde{S_i}}{\tilde{S_i}} = \mu_idt + \sum_{j=1}^{m}\sigma_{ij} dW_j \,\, , \,\, \forall i = 1,2,...,m$$
# 
# 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:
# 
# $$(\epsilon_1,\epsilon_2,..,\epsilon_m) = \begin{cases}
# (\epsilon_1 (\omega_1),\epsilon_2 (\omega_1),..,\epsilon_m (\omega_1)), \,\, \text{with probability } q_1, \\
# (\epsilon_1 (\omega_2),\epsilon_2 (\omega_2),..,\epsilon_m (\omega_2)), \,\, \text{with probability } q_2, \\
# (\epsilon_1 (\omega_3),\epsilon_2 (\omega_3),..,\epsilon_m (\omega_3)), \,\, \text{with probability } q_3 = (1-q_1-q_2)
# \end{cases}
# $$

# The equation describing the event tree (first equation of the notebook) becomes a discrete approximation of the above-stated stochastic process if we set:
# 
# $$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)$$
# 
# $$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)$$
# 
# $$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)$$
# 
# 
# 
# In the limit $T\to \infty$ (or alternatively $\Delta t \to 0$), the approximation approaches the continuous process.

# ### 3.) Investor's reservation purchase price for a European call option

# 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. 

# 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:

# $$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)]$$
# 
# 
# $$\text{subject to:   }\beta_0 \leq W_0 - C - \sum_{i=1}^{m} S_{i,0}[\alpha_{i,0} + |\alpha_{i,0}|],$$
# $$\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)|],$$
# $$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$$

# 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:
# 
# $$C^* = \text{maximize}_{({\mathbf{\alpha}_n(I_n),\beta_n(I_n))}_{I_n\in F_T,{n=1,2,..,T-1}}} \,\,\, C$$
# 
# $$\text{subject to: } E[U(W_T)] \geq U^*(W_0) $$
# 
# along with the other constraints mentioned in the optimization problem for the portfolio in which the option is purchased.

# ### 4.) Utility function

# 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:
# 
# $$U(W) = \frac{1-\gamma}{\gamma}\bigg(\frac{aW}{1-\gamma} + b\bigg)^\gamma \,;\, a>0$$
# 
# where $W > ((\gamma - 1)b)/a$. Note that the exponential and logarithmic utility functions are also members of the HARA-class.

# Another thing to consider is the Arrow-Pratt measure of the Absolute risk aversion. For the HARA-class, it is:
# 
# $$ARA(W) = \bigg(\frac{W}{1-\gamma} + \frac{b}{a}\bigg)^{-1}$$

# # Fusion Implementation

# 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.

# The various constants in the initialization of the Tree class are as follows:
# <b> 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 </b>.

# ### HARA utility as a Power cone:

# The objective function is:
# 
# $$\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$$
# 
# where $P(I_T)$ is the probability of a given path. We can re-write this as:
# 
# 
# $$\text{maximize   } \sum_{I_T\in F_T}P(I_T) \bigg(\frac{1-\gamma}{\gamma}\bigg) h(I_T)$$
# 
# $$\text{subject to :   } h(I_T)\leq \bigg(\frac{aW(I_T)}{1-\gamma} + b\bigg)^\gamma$$
# 
# 
# Now, the constraint can be expressed as a Power cone. However, there are a few possible cases:
# 
# 1.) $\gamma > 1$: In this case, the conic representation is:
# 
# $$\Bigg(h(I_T),1,\bigg(\frac{aW(I_T)}{1-\gamma} + b\bigg) \Bigg) \in \mathcal{P}_3^{1/\gamma}$$
# 
# 2.) $0 < \gamma < 1$:
# 
# $$\Bigg(\bigg(\frac{aW(I_T)}{1-\gamma} + b\bigg),1, h(I_T) \Bigg) \in \mathcal{P}_3^{\gamma}$$
# 
# 3.) $\gamma < 0$:
# 
# $$\Bigg(h(I_T),\bigg(\frac{aW(I_T)}{1-\gamma} + b\bigg),1  \Bigg) \in \mathcal{P}_3^{1/(1-\gamma)}$$

# ### Exponential utility as an exponential cone:

# The objective function in the case of exponential utility will be given by:
# 
# $$\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)
# $$
# 
# if the Absolute Risk aversion is nonzero. However, if ARA$(W) = 0$, then the objective function is simply:
# 
# $$\text{maximize   }E[U(W_T)] = \sum_{I_T\in F_T}P(I_T) W(I_T)$$
# 
# For the non-zero ARA$(W)$, we express the objective function as follows:
# 
# $$\text{maximize   } \sum_{I_T\in F_T}P(I_T) h(I_T)$$
# 
# $$\text{subject to:  } (1-h(I_T))\geq e^{-[\text{ARA}(W)]W(I_T)}$$
# 
# and the constraint is readily expressed as an exponential cone:
# 
# $$\big((1-h(I_T)),1,-[\text{ARA}(W)]W(I_T)\big) \in K_{exp}$$

# In[1]:


import numpy as np
from mosek.fusion import *
import sys
import matplotlib.pyplot as plt
import time


# In[2]:


class Tree:    
    #Instantiate the Tree class.
    def __init__(self,M,T,W_0,S,theta,S_v_0,r,U):
        self.M = M
        self.W_0 = W_0
        self.T = T
        self.S_v = S
        self.r = r    
        #Shape of the cost scaling matrix will give us the split index and the number of risky securities.
        self.s,self.n = S.shape
        self.cost = np.asarray([theta])
        self.S_v_0 = S_v_0
        #Parameters for the utility function.
        self.U = U
        self.S_v_t = S_v_0
        #This is to decide which utility will be used.
        self.util_dispatch = {'HARA': self.HARA_util,'EXP': self.exp_util}
        self.levels = []
    
    
    #Method to make all the levels, or the time steps in the tree.
    def level_make(self):
        #Number of risky securities, i.e. alpha{i,t=0} [Note that self.n is the number of risky securities]
        self.a0 = self.M.variable('a_0',[1,self.n])
        #Bonds at t=0.
        self.b0 = self.M.variable('b_0',1)
        #Making the first budget constraint.
        self.budg_ex = Expr.sub(Expr.sub(self.W_0,self.b0),
                                Expr.dot(self.S_v_0,Expr.mulElm(1+self.cost,self.a0)))
        self.root_budget = self.M.constraint('Budget_0',self.budg_ex,Domain.greaterThan(0.0))
        #This list will hold the level objects associated to this tree. (Please look at the sub-class: Level)
        a,b = self.a0,self.b0
        for i in range(1,self.T+1):
            #Appending level n, based on the level (n-1).
            self.levels.append(Level(i,self,a,b))
            a = self.levels[i-1].a_new
            b = self.levels[i-1].b_new
        #a_T is zero (). 

    
    #Method to make the HARA utility function and the corresponding Power-cone constraint.
    def HARA_util(self,W):
        self.h = self.M.variable('h',self.s**self.T,Domain.greaterThan(0.0))
        
        #HARA utility is only defined for when W > (gamma - 1)a/b.
        self.M.constraint(W,Domain.greaterThan(self.U[1]*(self.U[2]-1)/self.U[0]))

        I = Expr.constTerm([self.s**self.T],1.0)
        self.E1 = Expr.add(Expr.mul(self.U[0]/(1-self.U[2]),W),Expr.constTerm(W.getShape(),self.U[1]))
        
        #There are multiple cases for different values of gamma, modelled as follows
        #For gamma > 1:
        if self.U[2]>1.0:
            self.M.constraint('Obj_terms_HARA',Expr.hstack(self.h,I,self.E1),Domain.inPPowerCone(1/self.U[2]))
        else:
            if self.U[2]<0.0: #For gamma < 0
                self.M.constraint('Obj_terms_HARA',Expr.hstack(self.h,self.E1,I),Domain.inPPowerCone(1/(1-self.U[2])))
            else:             #For 0 < gamma < 1
                self.M.constraint('Obj_terms_HARA',Expr.hstack(self.E1,I,self.h),Domain.inPPowerCone(self.U[2]))    
        return(Expr.mul(((1-self.U[2])/self.U[2]),self.h))
    
    
    #Method to make the Exponential utility constraints.
    def exp_util(self,W):
        #Zero ARA(W) case
        if self.U == 0:
            return(W)
        #Nonzero ARA(W) case
        else:
            self.h = self.M.variable('h',self.s**self.T)
            I = Expr.constTerm([self.s**self.T],1.0)
            E_exp = Expr.mul(-self.U,W)
            self.M.constraint('Obj_terms_exp',Expr.hstack(Expr.sub(I,Expr.mul(self.U,self.h)),I,E_exp),
                              Domain.inPExpCone())
            return(self.h)


# The <i>Tree</i> class defined above has the sub-class <i>Level</i>. 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):
# 
# $$
# \begin{pmatrix}
# a\begin{bmatrix}
# S_{(1,2)} \\
# S_{(2,2)} \\
# S_{(3,2)} \\
# \end{bmatrix}\\
# b\begin{bmatrix}
# S_{(1,2)} \\
# S_{(2,2)} \\
# S_{(3,2)} \\
# \end{bmatrix}\\
# c\begin{bmatrix}
# S_{(1,2)} \\
# S_{(2,2)} \\
# S_{(3,2)} \\
# \end{bmatrix}\\
# \end{pmatrix}
# :
# \begin{pmatrix}
# \begin{bmatrix}
# \mathbf{\alpha_{(1,2)}},\beta_{(1,2)} \\
# \mathbf{\alpha_{(2,2)}},\beta_{(2,2)} \\
# \mathbf{\alpha_{(3,2)}},\beta_{(3,2)} \\
# \end{bmatrix} \\
# \begin{bmatrix}
# \mathbf{\alpha_{(1,2)}},\beta_{(1,2)} \\
# \mathbf{\alpha_{(2,2)}},\beta_{(2,2)} \\
# \mathbf{\alpha_{(3,2)}},\beta_{(3,2)} \\
# \end{bmatrix} \\
# \begin{bmatrix}
# \mathbf{\alpha_{(1,2)}},\beta_{(1,2)} \\
# \mathbf{\alpha_{(2,2)}},\beta_{(2,2)} \\
# \mathbf{\alpha_{(3,2)}},\beta_{(3,2)} \\
# \end{bmatrix} \\
# \end{pmatrix}  \longrightarrow
# \begin{pmatrix}
# \mathbf{\alpha_{(1,3)}},\beta_{(1,3)} \\
# \mathbf{\alpha_{(2,3)}},\beta_{(2,3)} \\
# \mathbf{\alpha_{(3,3)}},\beta_{(3,3)} \\
# \mathbf{\alpha_{(4,3)}},\beta_{(4,3)} \\
# \mathbf{\alpha_{(5,3)}},\beta_{(5,3)} \\
# \mathbf{\alpha_{(6,3)}},\beta_{(6,3)} \\
# \mathbf{\alpha_{(7,3)}},\beta_{(7,3)} \\
# \mathbf{\alpha_{(8,3)}},\beta_{(8,3)} \\
# \mathbf{\alpha_{(9,3)}},\beta_{(9,3)} \\
# \end{pmatrix}
# $$
# 
# Here, the first column vector is a <i>"vstack"</i> 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$.
# 
# The above-stated representation shows that we can <i>"repeat"</i> the variables of the previous level $s$ number of times and then it becomes quite easy to implement the budget constraints in <i>Fusion</i>. One can also extend the above shown method to a higher split index.

# In[3]:


#Each level corresponds to a time step in the event tree. This is a subclass of the 'Tree' class.
class Level(Tree):
    # l corresponds to the time step; a_old, b_old belong to (l-1) in the tree. 
    def __init__(self,l,Tree,a_old,b_old):
        if l==Tree.T :
            #If the current level is the final time period, then all risky securities are considered sold.
            self.a_new = Expr.constTerm([Tree.s**l,Tree.n],0.0)
            #Moreover, the final step's bonds are the wealth W(I_T); later used in the utility function.
            self.b_new = Tree.M.variable('W_T',Tree.s**l)
        else:
            #Risky securities for time step l.
            self.a_new = Tree.M.variable('a_{}'.format(l),[Tree.s**l,Tree.n])
            #Bonds in dollars for time step l.
            self.b_new = Tree.M.variable('b_{}'.format(l),Tree.s**l)    
        #Variable for the quadratic cone to implement the absolute value constraint.
        self.t_new = Tree.M.variable('t_{}'.format(l),[Tree.s**l,Tree.n],Domain.greaterThan(0.0))
        Tree.cost = np.repeat(Tree.cost,Tree.s,axis=0)

        #Repeating/Stacking the (l-1) level variables to implement the linear budget constraints. 
        A = Expr.repeat(a_old,Tree.s,0)
        B = Expr.repeat(b_old,Tree.s,0)
        #The price vector of the previous level, multiplied by the scale factors and then stacked vertically.
        Tree.S_v_t = np.vstack([np.multiply(j,Tree.S_v_t) for j in Tree.S_v])

        #Expressions involved in the budget constraints.
        self.bond_sub = Expr.sub(Expr.mul(B,(1+Tree.r)),self.b_new)
        self.secu_sub = Expr.sub(self.a_new,A)
        self.transact = Expr.add(self.secu_sub,Expr.mulElm(Tree.cost,self.t_new))
        self.secu_exp = Expr.mulDiag(self.transact,Tree.S_v_t.transpose())
        
        #The linear budget constraints.
        self.budg_constr = Tree.M.constraint('Budget_{}'.format(l),Expr.sub(self.bond_sub,self.secu_exp),
        Domain.greaterThan(0.0))
        #Quadratic cone for the absolute value requirement in the budget constraints.
        Tree.M.constraint('a{}_abs'.format(l),Expr.stack(2,self.t_new,self.secu_sub),Domain.inQCone())


# ### Now that we have the basic structure of the tree ready, we can easily make the <i>Fusion</i> model as follows...

# 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.

# In[4]:


#This function converts the path number from base 10 to base split-index. 
def base_rep(b,i,size):
    n = np.zeros(size)
    for j in range(size):
        n[j] = i%b
        if i//b == 0:
            break
        else:
            i = i//b
    return(np.flip(n))

#This function converts the path number to path id, for all paths.
def path_id_make(split_index,T):
    path_id = []
    for i in range(split_index**T):
        path_id.append(base_rep(split_index,i,T))
    return(np.asarray(path_id).astype(int))

#This function calculates: path probabilities as well as the price of the underlying stock.
def path_route_calc(path_id,split_index,q,*args):
    s = np.zeros(split_index)
    for p in path_id:
        s[p] = s[p] + 1
    if args:
        path_S1 = []
        for j in range(split_index):
            path_S1.append((args[0][j])**s[j])
        path_S1T = np.prod(np.asarray(path_S1))
        #Returning price of the stock.
        return(path_S1T)
    else:
        path_prob = np.prod(q**s)
        #Returning probability of the path.
        return(path_prob)


#These functions will be used for the calculation of the price-scaling factors (i.e. a,b and c).
#They follow the discrete approximation of the continuous price vector process.
def price_vector_z(sigma,eps,dt,q):
    E = np.exp(np.matmul(eps,sigma)*np.sqrt(dt))
    return(np.matmul(q,E))

def price_vector_abc(sigma,eps,dt,q,mu,Z):
    E1 = np.multiply(np.exp(np.matmul(eps,sigma)*np.sqrt(dt)),1/Z)
    E2 = np.repeat([np.exp(mu*dt)],eps.shape[0],axis = 0)
    return(np.multiply(E2,E1))


# ### Portfolio without call option:

# In[5]:


def Portfolio(port_params,util_type = 'HARA',wrtLog = True):
    T,W_0,S,theta,S_v_0,r,U,q = port_params
    M = Model('PORTFOLIO')
    if wrtLog:
        M.setLogHandler(sys.stdout)
    #Instantiating the Tree class.    
    Tree_1 = Tree(M,T,W_0,S,theta,S_v_0,r,U)
    #Making all the levels.
    Tree_1.level_make()
    #Making the conic constraints for the utility function.
    H = Tree_1.util_dispatch[util_type](Tree_1.levels[T-1].b_new)
    #Making the path ids (see function def).
    path_ids = path_id_make(Tree_1.s,T)
    #Calculating the path probabilities (see function def).
    path_probs = [path_route_calc(p,Tree_1.s,q) for p in path_ids]
    #Objective function
    Obj = Expr.dot(H,path_probs)
    M.objective('PORTFOLIO_OBJ',ObjectiveSense.Maximize,Obj)
    M.solve()
    #The maximal expected utility
    utility_W0 = M.primalObjValue()
    ut_time = M.getSolverDoubleInfo('optimizerTime')
    ut_iter = M.getSolverIntInfo('intpntIter')
    return(M,Tree_1,utility_W0,path_ids,Obj,[ut_iter,ut_time])


# ### Portfolio with call option:

# 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).

# In[6]:


def Option_Portfolio(port_params,K,util_type = 'HARA',writeLog = True,solver_info=False):
    T,W_0,S,theta,S_v_0,r,U,q = port_params
    #We start by solving the problem for the no option case, so that we have U*(W_0).
    M,Tree_1,utility_W0,path_ids,Obj,obj_info = Portfolio(port_params,util_type = util_type,wrtLog = writeLog)
    #This is the price of the underlying stock at the Time horizon (for each path, i.e. s**T paths)
    path_Svs = np.asarray([path_route_calc(p,S.shape[0],q,S[:,0]) for p in path_ids])
    #max((Stock_price(I_T) - K),0)
    call_profit = ((path_Svs - K) + abs(path_Svs - K))/2    
    #Adding a new variable to the model: the reservation purchase price of the option.
    #Note: the call option is bought only on one underlying security.
    Call = M.variable('Call_Price',1,Domain.inRange(0.0,S_v_0[0]))
    #Updating the root constraint, making W_0 into (W_0 - C).
    Tree_1.root_budget.update(Expr.neg(Call),Call)
    #Updating the budget constraints in the final time-period by adding max((Stock_price(I_T) - K),0).
    Tree_1.levels[T-1].budg_constr.update(call_profit.tolist())
    #Now, we add the constraint E[U(W_T)] >= U*(W_0), to our model.
    M.constraint('E_geq_U',Expr.sub(Obj,utility_W0),Domain.greaterThan(0.0))
    #Objective function: given by the reservation purchase price itself
    M.objective('Call_price_OBJECTIVE',ObjectiveSense.Maximize,Call)
    M.solve()
    if solver_info:
        price_time = M.getSolverDoubleInfo('optimizerTime')
        price_iter = M.getSolverIntInfo('intpntIter')
        n_cons = M.getSolverIntInfo('optNumcon')
        n_vars = M.getSolverIntInfo('optNumvar')
        return(M.primalObjValue(),n_cons,n_vars,obj_info[0],obj_info[1],price_iter,price_time)
    else:
        return(M.primalObjValue())


# Now, to demonstrate the model, we solve it for two simple cases. The values taken below are the same as mentioned in the paper (<a href="https://www.sciencedirect.com/science/article/pii/S0168927498001044">Andersen et. al.</a>). 

# In[7]:


T = 5                                                              #Number of periods
sigma = np.asarray([0.2])                                          #Sigma matrix
q = np.ones(3)/3                                                   #Probabilities
eps = np.asarray([[np.sqrt(2)],[-1/np.sqrt(2)],[-1/np.sqrt(2)]])   #Eps matrix (see text)
T_horizon = 1/4                                                    #Time horizon (in years)
dT = T_horizon/T                                                   
theta = [0.005]                                                    #Transaction costs
S_v_0 = [1.0]                                                      #Initial price vector
r = (1.06**dT) - 1                                                 #Constant interest rate
W_0 = 1.0                                                          #Initial wealth
gamma = 0.3                                                        #Gamma parameter for HARA utility
b = 1                                                              #b parameter for HARA utility (see text)
c = 0.2                                                            #Absolute risk aversion ARA(W)
a = b/((1/c) +  (W_0/(gamma-1)))
K = 1                                                              #Strike Price.
mu = 0.15                                                          #Drift of the stock price

Z = price_vector_z(sigma,eps,dT,q)
S_coeffs = np.expand_dims(price_vector_abc(sigma,eps,dT,q,mu,Z),axis=1)
    
util_paras = np.asarray([a,b,gamma])
input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,util_paras,q]

#Using HARA utility
print('\n\nCall option price = {}'.format(Option_Portfolio(input_pars,K)))


# In[8]:


#Using Exponential utility (everything else is the same as above)
input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]
print('\n\nCall option price = {}'.format(Option_Portfolio(input_pars,K,util_type='EXP')))


# # Test-cases:

# ### Reservation price sensitivity for the choice of utility function:

# <b>(a.) Dependency of the reservation purchase price on $\gamma$, $ARA(W_0)$ and $\overline{T}$.</b> (Consider Table 3, <a href="https://www.sciencedirect.com/science/article/pii/S0168927498001044">Andersen et. al.</a>)

# In[9]:


gamma = [-4.0,-2.0,-0.9,-0.3,0.3,0.6]                              #Gamma parameter for HARA utility

print('ARA = 0.2 ; Time Horizon = 1/4 year')
print('{0:^6}  {1:^9}'.format('gamma','C'))
for g in gamma:
    a = b/((1/c) +  (W_0/(g-1)))
    util_paras = np.asarray([a,b,g])
    input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,util_paras,q]
    #Using HARA utility
    print('{0:^ 5.1f}  {1:^9.7f}'.format(g,Option_Portfolio(input_pars,K,writeLog=False)))
input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]
#Using Exponential utility
print('{0:^5}  {1:^9.7f}'.format('exp',Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))


# In[10]:


c = 1.0                                                            #Absolute risk aversion ARA(W)

print('ARA = 1.0 ; Time Horizon = 1/4 year')
print('{0:^6}  {1:^9}'.format('gamma','C'))
for g in gamma[0:4]:
    a = b/((1/c) +  (W_0/(g-1)))
    util_paras = np.asarray([a,b,g])
    input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,util_paras,q]
    #Using HARA utility
    print('{0:^ 5.1f}  {1:^9.7f}'.format(g,Option_Portfolio(input_pars,K,writeLog=False)))
input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]
#Using Exponential utility
print('{0:^5}  {1:^9.7f}'.format('exp',Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))


# In[11]:


T_horizon = 9                                                      #Time horizon (in years)
dT = T_horizon/T                             
Z = price_vector_z(sigma,eps,dT,q)
S_coeffs = np.expand_dims(price_vector_abc(sigma,eps,dT,q,mu,Z),axis=1)    
r = (1.06**dT) - 1                                                 #Constant interest rate
c = 0.2                                                            #Absolute risk aversion ARA(W)

print('ARA = 0.2 ; Time Horizon = 9 years')
print('{0:^6}  {1:^9}'.format('gamma','C'))
for g in gamma:
    a = b/((1/c) +  (W_0/(g-1)))
    util_paras = np.asarray([a,b,g])
    input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,util_paras,q]
    #Using HARA utility
    print('{0:^ 5.1f}  {1:^9.7f}'.format(g,Option_Portfolio(input_pars,K,writeLog=False)))
input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]
#Using Exponential utility
print('{0:^5}  {1:^9.7f}'.format('exp',Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))


# In[12]:


c = 1.0                                                            #Absolute risk aversion ARA(W)

print('ARA = 1.0 ; Time Horizon = 9 years')
print('{0:^6}  {1:^9}'.format('gamma','C'))
for g in gamma[0:4]:
    a = b/((1/c) +  (W_0/(g-1)))
    util_paras = np.asarray([a,b,g])
    input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,util_paras,q]
    #Using HARA utility
    print('{0:^ 5.1f}  {1:^9.7f}'.format(g,Option_Portfolio(input_pars,K,writeLog=False)))
input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]
#Using Exponential utility
print('{0:^5}  {1:^9.7f}'.format('exp',Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))


# <b>(b.) Dependency of the reservation purchase price on the initial level of Relative Risk Aversion, $RRA(W_0)$</b> (Consider Table 4, <a href="https://www.sciencedirect.com/science/article/pii/S0168927498001044">Andersen et. al.</a>)

# In[13]:


T = 5                                                              #Number of periods
T_horizon = 1/4                                                    #Time horizon (in years)
dT = T_horizon/T                                                   
r = (1.06**dT) - 1                                                 #Constant interest rate
W_0 = [1.0,4.0,8.0]                                                #Initial wealth
gamma = -1.0                                                       #Gamma parameter for HARA utility
b = 1                                                              #b parameter for HARA utility (see text)
c = 0.2                                                            #Absolute risk aversion ARA(W)
Z = price_vector_z(sigma,eps,dT,q)
S_coeffs = np.expand_dims(price_vector_abc(sigma,eps,dT,q,mu,Z),axis=1)

print('ARA = 0.2 ; Time Horizon = 1/4 year')
print('{0:^5}  {1:^5}  {2:^9}'.format('W_0','RRA','C'))
for w0 in W_0:
    a = b/((1/c) +  (w0/(gamma-1)))
    util_paras = np.asarray([a,b,gamma])
    input_pars = [T,w0,S_coeffs,theta,S_v_0,r,util_paras,q]
    #Using HARA utility
    print('{0:^5.0f}  {1:^5.1f}  {2:^9.7f}'.format(w0,c*w0,Option_Portfolio(input_pars,K,writeLog=False)))
input_pars = [T,1.0,S_coeffs,theta,S_v_0,r,c,q]
#Using Exponential utility
print('{0:^5}  {1:^5}  {2:^9.7f}'.format('exp','',Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))


# In[14]:


T_horizon = 9                                                      #Time horizon (in years)
dT = T_horizon/T
r = (1.06**dT) - 1                                                 #Constant interest rate
Z = price_vector_z(sigma,eps,dT,q)
S_coeffs = np.expand_dims(price_vector_abc(sigma,eps,dT,q,mu,Z),axis=1)

print('ARA = 0.2 ; Time Horizon = 9 years')
print('{0:^5}  {1:^5}  {2:^9}'.format('W_0','RRA','C'))
for w0 in W_0:
    a = b/((1/c) +  (w0/(gamma-1)))
    util_paras = np.asarray([a,b,gamma])
    input_pars = [T,w0,S_coeffs,theta,S_v_0,r,util_paras,q]
    #Using HARA utility
    print('{0:^5.0f}  {1:^5.1f}  {2:^9.7f}'.format(w0,c*w0,Option_Portfolio(input_pars,K,writeLog=False)))
input_pars = [T,1.0,S_coeffs,theta,S_v_0,r,c,q]

print('{0:^5}  {1:^5}  {2:^9.7f}'.format('exp','',Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))


# <b>(c.) Convergence of the reservation purchase price.</b> (Consider Table 5, <a href="https://www.sciencedirect.com/science/article/pii/S0168927498001044">Andersen et. al.</a>)

# In[15]:


T_horizon = 1/4                                                    #Time horizon (in years)
W_0 = 1.0                                                          #Initial wealth
gamma = 0.3                                                        #Gamma parameter for HARA utility
b = 0.0                                                            #b parameter for HARA utility (see text)
c = 0.7                                                            #Absolute risk aversion ARA(W)
a = 1 - gamma
util_paras = np.asarray([a,b,gamma])

C_trio = []
print('{0:^3}  {1:12}  {2:12}  {3:^12}'.format(' ','-exp(-0.7W)','(W^0.3)/0.3',' '))
print('{0:^3}  {1:^12}  {2:^12}  {3:^12}'.format('T','C_exp','C_pow','CRR'))
for T in np.arange(1,12):
    dT = T_horizon/T
    r = (1.06**dT) - 1
    Z = price_vector_z(sigma,eps,dT,q)
    S_coeffs = np.expand_dims(price_vector_abc(sigma,eps,dT,q,mu,Z),axis=1)
    #C_pow : Power utility function.
    input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,util_paras,q]
    C_pow = Option_Portfolio(input_pars,K,writeLog=False)
    #CRR : Call price in a frictionless market 
    input_pars = [T,W_0,S_coeffs,[0.0],S_v_0,r,util_paras,q]
    CRR = Option_Portfolio(input_pars,K,writeLog=False)
    #C_exp : Exponential utility function.
    input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]
    C_exp = Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)
    print('{0:^3d}  {1:^12.10f}  {2:^12.10f}  {3:^12.10f}'.format(T,C_exp,C_pow,CRR))
    C_trio.append([C_exp,C_pow,CRR])


# <b>(d.) Reservation purchase price dependence on absolute risk aversion (Exponential utility).</b> (Consider Table 6, <a href="https://www.sciencedirect.com/science/article/pii/S0168927498001044">Andersen et. al.</a>)

# In[16]:


T = 9                                                              #Number of periods
dT = T_horizon/T                                                   
theta = [0.005]                                                    #Transaction costs
r = (1.06**dT) - 1                                                 #Constant interest rate
c = [0.1,0.2,0.4,0.8,1.0,2.0,4.0]                                  #Absolute risk aversion ARA(W)
Z = price_vector_z(sigma,eps,dT,q)
S_coeffs = np.expand_dims(price_vector_abc(sigma,eps,dT,q,mu,Z),axis=1)    

print('{0:^8}  {1:^12}'.format('ARA(W_0)','C[ARA(W_0)]'))
C_ara = []
for kappa in c:
    input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,kappa,q]
    C_exp = Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)
    print('{0:^8.2f}  {1:^12.10f}'.format(kappa,C_exp))
    C_ara.append(C_exp)


# <b>(e.) Reservation purchase price dependence on absolute risk aversion (Power utility).</b> (Consider Table 7, <a href="https://www.sciencedirect.com/science/article/pii/S0168927498001044">Andersen et. al.</a>)

# In[17]:


eta = np.linspace(0.1,0.9,9)                                        #Risk Aversion parameter
a = 1
b = 0.0

print('{0:^5}  {1:^12}'.format('eta','C[ARA(W_0)]'))
C_eta = []
for n in eta:
    gamma = 1-n
    util_paras = np.asarray([a,b,gamma])
    input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,util_paras,q]
    C_pow = Option_Portfolio(input_pars,K,writeLog=False)
    print('{0:^5.2f}  {1:^12.10f}'.format(n,C_pow))
    C_eta.append(C_pow)


# <b>(f.) Reservation purchase price dependence on initial wealth (Power utility with $\eta = 0.7$).</b> (Consider Table 8, <a href="https://www.sciencedirect.com/science/article/pii/S0168927498001044">Andersen et. al.</a>)

# In[18]:


eta = 0.7                                                            #Risk Aversion parameter
W_0 = [0.25,0.50,1.0,2.0,4.0,8.0]                                    #Initial wealth
gamma = 1-eta
util_paras = np.asarray([a,b,gamma])

print('{0:^5}  {1:^12}'.format('W_0','C_pow(W_0)'))
C_w0 = []
for w0 in W_0:
    input_pars = [T,w0,S_coeffs,theta,S_v_0,r,util_paras,q]
    C_pow = Option_Portfolio(input_pars,K,writeLog=False)
    print('{0:^5.2f}  {1:^12.10f}'.format(w0,C_pow))
    C_w0.append(C_pow)


# ### The effect of diversification opportunities on the reservation purchase price.

# 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 <a href="https://www.sciencedirect.com/science/article/pii/S0168927498001044">Andersen et. al.</a>. For this example, we have:
# 
# $$
# \sigma^{T} = \begin{pmatrix}
# \sigma_{11}\, ,\,\sigma_{21} \\
# \sigma_{12}\, ,\,\sigma_{22} \\
# \end{pmatrix} = \begin{pmatrix}
# \sigma_{11}\, ,\, 0.2\rho\\
# 0.0\, , \, \sqrt{(0.2)^2 - \sigma_{21}^2} \\
# \end{pmatrix}
# $$
# 
# where $\rho$ is the correlation between logarithms of the two securities. The $\epsilon$ matrix is given by, 
# 
# $$\epsilon = \begin{pmatrix}
# \epsilon_1 (\omega_1),\epsilon_2 (\omega_1) \\
# \epsilon_1 (\omega_2),\epsilon_2 (\omega_2) \\
# \epsilon_1 (\omega_3),\epsilon_2 (\omega_3)
# \end{pmatrix} = \begin{pmatrix}
# \sqrt{2}\,\,,\,\,0 \\
# -1/\sqrt{2},\sqrt{3/2} \\
# -1/\sqrt{2},\sqrt{3/2}
# \end{pmatrix}
# $$

# 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.
# 
# <b>Note</b>: 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.

# <b>(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.</b> (Consider Table 10, <a href="https://www.sciencedirect.com/science/article/pii/S0168927498001044">Andersen et. al.</a>)

# In[19]:


#No transaction costs (Friction-less market)
W_0 = 1.0                                                          #Initial wealth
gamma = 0.3                                                        #Gamma parameter for HARA utility
b = 1                                                              #b parameter for HARA utility (see text)
c = 0.2                                                            #Absolute risk aversion ARA(W)
a = b/((1/c) +  (W_0/(gamma-1)))
util_paras = np.asarray([a,b,gamma])
K = 1                                                              #Strike Price.

#Parameters for the 2 risky securities under consideration:
#The sigma matrix entries for the first security are same as before...
sig11 = 0.2                                                        #sigma_{1,1}
sig12 = 0.0                                                        #sigma_{1,2}
#Epsilon matrix (Note the extra column for the second risky security)
eps = np.asarray([[np.sqrt(2),0],[-1/np.sqrt(2),np.sqrt(3/2)],[-1/np.sqrt(2),-np.sqrt(3/2)]])
S_v_0 = [1.0,1.0]                                                  #Initial price vector
mu = np.asarray([0.15,0.15])                                       #Drift of the stock price
rho = [1,0.5,0.0,-0.5,-0.9]                                        #Correlation

theta = [0.0,0.0]                                                  #Transaction costs
print('theta_1 = {0:5.3f} ; theta_2 = {1:5.3f}'.format(theta[0],theta[1]))
print('{0:^5}  {1:^12}'.format('rho','C(S1,S2)'))
for p in rho:
    sig21 = p*0.2
    #Sigma matrix
    sigma = np.asarray([[sig11,sig21],[sig12,np.sqrt((0.2**2) - (sig21**2))]])
    Z = price_vector_z(sigma,eps,dT,q)
    S_coeffs = price_vector_abc(sigma,eps,dT,q,mu,Z)
    input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]
    print('{0:^ 4.1f}  {1:^12.10f}'.format(p,Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))


# In[20]:


theta = [0.005,0.0]                                                 #Transaction costs
print('theta_1 = {0:5.3f} ; theta_2 = {1:5.3f}'.format(theta[0],theta[1]))
print('{0:^5}  {1:^12}'.format('rho','C(S1,S2)'))
for p in rho:
    sig21 = p*0.2
    #Sigma matrix
    sigma = np.asarray([[sig11,sig21],[sig12,np.sqrt((0.2**2) - (sig21**2))]])
    Z = price_vector_z(sigma,eps,dT,q)
    S_coeffs = price_vector_abc(sigma,eps,dT,q,mu,Z)
    input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]
    print('{0:^ 4.2f}  {1:^12.10f}'.format(p,Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))


# In[21]:


theta = [0.005,0.005]                                               #Transaction costs
print('theta_1 = {0:5.3f} ; theta_2 = {1:5.3f}'.format(theta[0],theta[1]))
print('{0:^5}  {1:^12}'.format('rho','C(S1,S2)'))
for p in rho:
    sig21 = p*0.2
    #Sigma matrix
    sigma = np.asarray([[sig11,sig21],[sig12,np.sqrt((0.2**2) - (sig21**2))]])
    Z = price_vector_z(sigma,eps,dT,q)
    S_coeffs = price_vector_abc(sigma,eps,dT,q,mu,Z)
    input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]
    print('{0:^ 4.2f}  {1:^12.10f}'.format(p,Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))


# <b>(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$.</b> (Consider Table 11, <a href="https://www.sciencedirect.com/science/article/pii/S0168927498001044">Andersen et. al.</a>)

# In[ ]:


rho = -0.9                                                          #Correlation
sig21 = rho*0.2
sigma = np.asarray([[sig11,sig21],[sig12,np.sqrt((0.2**2) - (sig21**2))]])
Z = price_vector_z(sigma,eps,dT,q)
S_coeffs = price_vector_abc(sigma,eps,dT,q,mu,Z)

theta_list = [0.0008,0.0016,0.003,0.006,0.01,0.02,0.05,0.1]
print('{0:^5}  {1:^12}'.format('theta','C(S1,S2)'))
for theta_i in theta_list:
    #We take the transaction costs for both securities to be equal.
    theta = [theta_i,theta_i]
    input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]
    print('{0:^5.4f}  {1:^12.10f}'.format(theta_i,Option_Portfolio(input_pars,K,util_type='EXP',writeLog=False)))


# <b>(i.) Computational efficiency for different values of T (time periods).</b> (Consider Table 12, <a href="https://www.sciencedirect.com/science/article/pii/S0168927498001044">Andersen et. al.</a>)

# 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 <i>no-option</i> model as well as the <i>option-price</i> model), time. (time taken, in seconds for both models).

# In[ ]:


theta = [0.005,0.005]
print('{0:^2}  {1:^8}  {2:^10}  {3:^11}  {4:^11}'.format('','','','No-option','Option-price'))
print('{0:>2}  {1:>8}  {2:>10}  {3:^3}  {4:^6}  {5:^3}  {6:^6}'.
      format('T','cons.','vars.','it.','time','it.','time'))
total_info_arr = []
for T in range(1,10):
    try:
        T_horizon = 1/4
        dt = T_horizon/T
        r = (1.06**dT) - 1
        Z = price_vector_z(sigma,eps,dT,q)
        S_coeffs = price_vector_abc(sigma,eps,dT,q,mu,Z)
        input_pars = [T,W_0,S_coeffs,theta,S_v_0,r,c,q]
        call,n_cons,n_var,ut_it,ut_time,price_it,price_time = Option_Portfolio(input_pars,K,util_type='EXP',
                                                                               writeLog=False,solver_info=True)
        print('{0:>2d}  {1:>8d}  {2:>10d}  {3:>3d}  {4:>6.3f}  {5:>3d}  {6:>6.3f}'.
              format(T,n_cons,n_var,ut_it,ut_time,price_it,price_time))
        total_info_arr.append([call,n_cons,n_var,ut_it,ut_time,price_it,price_time])
    except MemoryError as e:
        print(e)
        break
    except SolutionError as s:
        print(s)
    except OptimizerError as e:
        print(e)
        break


# 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.

# <a rel="license" href="http://creativecommons.org/licenses/by/4.0/"><img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by/4.0/80x15.png" /></a><br />This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by/4.0/">Creative Commons Attribution 4.0 International License</a>. The **MOSEK** logo and name are trademarks of <a href="http://mosek.com">Mosek ApS</a>. 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).