% Applications of Linear Programming to Financial Mathematics
%
% Problem: find prices for hedging portfolios in incomplete markets
% Method: optimize the price of the portfolio given constraints expressed as
% linear inequalities that mean "cover all liabilities in all modeled states."
% Tool: Use the "glpk" package from the GNU software archive:
%    https://www.gnu.org/software/glpk/
%  Download, compile, and install the package on your computer.
%  Documentation on how to use it within Octave is available online:
%    https://octave.org/doc/v4.4.1/Linear-Programming.html
%
% DATE: 2020-12-17
% AUTHOR: M.V.Wickerhauser
% 
% From the trinomial example in Appendix B of Van der Hoek and Elliot, p.256
% Section B.2  Solutions to Incomplete Markets
% Define, by its payoff, the derivative to be hedged in an incomplete market:
K = 78                          % Strike price
payoff = @(s) max(s-K,0)        % European call option
%
% First, use constrained optimization to find a bid-ask interval.
%
% Parts (b), Superreplication and (c), Subreplication
% Define the optimization problem:
%  Find portfolio x such that c'*x is optimal given constraints on A*x versus b
S0 = 80;                        % Spot price of the underlying asset
c = [S0 1]'                     % Cost of the hedging portfolio
R = 1.01                        % Riskless return to expiry.
S1 = [85 80 75]'                % Modeled prices for the risky asset at expiry
b = payoff(S1)                  % Payoffs in the modeled states
A = [S1 R*ones(size(S1))]       % Constraint matrix
% Set some parameters for the linear programming function glpk()
vartype = "CC";                 % Variable type: x(j) is a Continuous variable, j=1,2
param.msglev=1;                 % Struct member for verbosity in glpk()
param.itlim=100;                % Struct member for maximum number of iterations to use
lb=[-S0 -S0]'                   % Lower bounds on x: sell short, or bank, no more than S0
ub=[S0 S0]'                     % Upper bounds on x: buy, or borrow, no more than S0
%%%  Use the function glpk() to compute the bid-ask spread
% (b) Superreplication, the derivative seller's perspective:
sellctype = "LLL";              % Constraint type: A*x(i) is Lower bounded by b(i), i=1,2,3
sellsense = 1;                  % Optimization direction for c'*x. "1" means "min"
[xmin, ask, status, extra] = glpk (c, A, b, lb, ub, sellctype, vartype, sellsense, param);
xmin                            % Cost-minimizing hedge portfolio
ask                             % "Ask" price  C^* for the Call: minimum cost to buy the hedge portfolio
% (c) Subreplication, the derivative buyer's perspective:
buyctype = "UUU";               % Constraint type: A*x(i) is Upper bounded by b(i), i=1,2,3
buysense = -1;                  % Optimization direction for c'*x: "-1" means "max"
[xmax, bid, status, extra] = glpk (c, A, b, lb, ub, buyctype, vartype, buysense, param);
xmax                            % Cost-maximizing hedge portfolio
bid                             % "Bid" price C_* for the Call: maximum cost before the hedge portfolio is better 
% Publish the bid-ask spread:
[bid, ask]
%
% Second, append another asset to get a complete market.
%
% Part (a) Complete the Market
% Find another derivative to complete the market
C0 = 2                          % Spot cost of the extra derivative
K0 = 80                         % Strike price of the extra derivative, must differ from K
extra = @(s) max(s-K0,0)        % Payoff of the extra derivative
Aa = [A extra(S1)]              % Augmented constraint matrix Aa is now invertible
ca = [c; C0]                    % Augmented spot cost of the hedging portfolio
% Solve for the unique hedging portfolio and find its price
hedge = Aa\b                    % Augmented hedge portfolio in the completed market
price = ca'*hedge               % Price of the hedge portfolio including the extra derivative


% Example from Chapter 8

K=88; payoff=@(s)max(s-K,0); % European-style Call at strike K
S0=90; % Spot price of the risky asset
c=[1;S0]; % q: spot prices of all market assets
R=1.02; % B(T): riskless return to expiry.
S1=[85;90;95]; % S(T): modeled payoffs for S
b=payoff(S1);  % C(T): modeled payoffs for Call
A=[R*ones(size(S1)),S1]; % A: market transposed into GLPK format

vartype="CC";    % ==> x(j) is Continuous, j=1,2
param.msglev=1;  % ==> use a low verbosity level
huge=1000;
param.itlim=huge; % ==> huge numberof iterations allowed
ub=[huge*S0;huge]; % huge Upper bounds on x
lb= -ub;           % huge Lower bounds on x

sellctype="LLL"; % Lower constraint type A*x(j)>=b(j), j=1,2,3
sellsense=1; % Optimization direction for c'*x: "1" ==> "min"
[hs,ask]=glpk(c,A,b,lb,ub,sellctype,vartype,sellsense,param);
hs % = [ -58.33333; 0.70000]; Cost-minimizing hedge portfolio 
ask % = 4.6667; minimum cost to superreplicate the Call

buyctype="UUU"; % Upper constraint type A*x(j)=<b(j), j=1,2,3
buysense=-1; % Optimization direction for c'*x: "-1" ==> "max"
[hb,bid]=glpk(c,A,b,lb,ub,buyctype,vartype,buysense,param);
hb % = [ -86.2745; 1.0000]; Cost-maximizing hedge portfolio
bid % = 3.7255; maximum cost to subreplicate the Call


% Exercise from Chapter 8

a=[10,15,20,25,30];  a0=20;
b=[17,15,12,10,7];   b0=12;
R=1.02;
A=[R R R R R; a; b];
q=[1; a0; b0];
A,q
rref([A R*q])
A\q*R
A\q

sellctype="LLLLL"; % Lower constraint type A'*x(j)>=bb(j), j=1:5
sellsense=1; % Optimization direction for q'*x: "1" ==> "min"
buyctype="UUUUU"; % Upper constraint type A'*x(j)=<bb(j), j=1:5
buysense=-1; % Optimization direction for q'*x: "-1" ==> "max"
vartype="CCC";   % ==> x(j) is Continuous, j=1:3
param.msglev=1;  % ==> use a low verbosity level
huge=1000; infty=huge*(abs(a0)+abs(b0)); % huge and huger
param.itlim=huge; % ==> huge maximum number of iterations
lb=[-infty; -huge; -huge]; % huge Lower bounds on x
ub=[ infty;  huge;  huge]; % huge Upper bounds on x

Ka=20; ba=max(a-Ka,0); % C(T): payoff for "a" Call
[hs,ask]=glpk(q,A',ba,lb,ub,sellctype,vartype,sellsense,param)
[hb,bid]=glpk(q,A',ba,lb,ub,buyctype,vartype,buysense,param)

Kb=13; bb=max(Kb-b,0); % P(T): payoff for "b" Put
[hs,ask]=glpk(q,A',bb,lb,ub,sellctype,vartype,sellsense,param)
[hb,bid]=glpk(q,A',bb,lb,ub,buyctype,vartype,buysense,param)
