-
Notifications
You must be signed in to change notification settings - Fork 5
/
finsolve.m
218 lines (202 loc) · 5.82 KB
/
finsolve.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
% FINSOLVE Solves continuous time asset pricing problems
% USAGE
% [c,V]=finsolve(model,fspace,alg,s,N);
% INPUTS
% model : a model structure (see below)
% fspace : name of projection space structure
% alg : algorithm used (lines, explicit, implicit, CN, or stiff)
% s : collocation nodal values of the state
% N : number of time steps
% OUTPUTS
% c : value function projection coefficients n x 1
% (or n x N+1 with keepall option)
% V : asset value n x 1
% (or n x N+1 with keepall option)
%
% The model structure has the following fields:
% model.func - the name of the model function file
% model.T - the time to maturity of the asset
% model.american : 1 if there is an early exercise feature for buyer
% -1 if there is an early exercise feature for seller
% model.barrier - 2x(2+d) matrix: rows specify lower and upper barriers
% model.params - additional parameters to be passed to model.func
% To price barrier options model.barrier should have the form
% [R_a alpha_a beta_a;
% R_b alpha_b beta_b];
% where R_a is paid anytime beta_a*S<=alpha_a and
% R_b is paid anytime beta_b*S<=alpha_b.
%
% The model function file must have the following syntax:
% function out=func(flag,S,additional parameters);
% switch flag
% case 'rho'
% out= instantaneous risk-free interest rate
% case 'mu'
% out= drift on the state process
% case 'sigma'
% out = volatility on the state process
% case 'delta'
% out = the payout flow (dividend) on the derivative asset
% case 'V0'
% out = exercise value of the asset
% end
%
% Note: The text documentation has an argument, t, which is not
% incorporated in the current version and is not used in the examples.
%
% USER OPTIONS (SET WITH OPSET)
% keepall : keeps coefficient values (c) for all time steps
% (this may be limited by available memory)
% payatT : for barrier options, 1 if rebate paid at expiration
% Revision 9/1/2010 replaced the use of lusolve with the pair luget/lusol
% Copyright (c) 1997-2010, Paul L. Fackler & Mario J. Miranda
% paul_fackler@ncsu.edu, miranda.4@osu.edu
function [c,V]=finsolve(model,fspace,alg,s,N)
%clc
keepall = optget('finsolve','keepall',0);
payatT = optget('finsolve','payatT',0);
if ~exist('alg','var') || isempty(alg), alg = 'lines'; end
if ~exist('s','var') || isempty(s), s = funnode(fspace); end
if ~exist('N','var') || isempty(N), N = 1; end
mffile = model.func;
T = model.T;
params = model.params;
n = prod(fspace.n);
d = fspace.d;
S = gridmake(s); % expand the grid
if isfield(model,'american'), american=model.american;
else american=0;
end
% Set up barier problem, if appropriate
if isfield(model,'barrier')
barrier=1;
temp=model.barrier;
Ra=temp(1,1); Rb=temp(2,1); % barrier rewards
Binda=S*temp(1,3:end)'<=temp(1,2); % points less than lower barrier
Bindb=S*temp(2,3:end)'>=temp(2,2); % points greater than upper barrier
else
barrier=0;
end
% Compute collocation matrix
bases=funbasx(fspace,s,[0;1;2]);
B=ctbasemake(feval(mffile,'mu',S,params{:}),bases,1);
B=B+ctbasemake(feval(mffile,'sigma',S,params{:}),bases,2);
rho=ctbasemake(feval(mffile,'rho',S,params{:}),bases,0);
B=B-rho;
if barrier
if payatT
B(Binda,:)=-rho(Binda,:);
B(Bindb,:)=-rho(Bindb,:);
else
B(Binda,:)=0;
B(Bindb,:)=0;
end
end
clear rho
Phi=bases.vals(1,d:-1:1);
bases.vals([2;3],:)=[];
% Get dividend information
b=feval(mffile,'delta',S,params{:});
hasdiv=~isempty(b);
% Get initial (terminal) value and compute initial coefficients
V0=feval(mffile,'V0',S,params{:});
if barrier
V0(Binda)=Ra;
V0(Bindb)=Rb;
end
c0=ckronxi(Phi,V0);
% Define matrix operators (or call ODE solver for stiff method)
Delta=T/N;
mass=1;
tensor=0;
switch alg
case 'lines'
B=ckronxi(Phi,B); % inv(Phi)*B
A=expm(full(Delta*B));
if hasdiv, a=(A-speye(n))*(B\ckronxi(Phi,b)); end
mass=0;
case 'explicit'
A=ckron(Phi)+Delta*B;
if hasdiv, a=Delta*b;end
tensor=1;
mass=0;
case 'implicit'
A=ckron(Phi);
M=A-Delta*B;
if hasdiv, a=Delta*b;end
MLU=luget(M);
case {'CN','trapezoid'}
B=(Delta/2)*B;
A=ckron(Phi);
M=A-B;
A=A+B;
if hasdiv, a=Delta*b; end
MLU=luget(M);
case {'stiff','ode'}
if alg(1)=='o'
solver=@ode45;
else
solver=@ode15s;
end
Phi=ckron(Phi);
if ~hasdiv, b=[]; end
options=odeset('Jacobian','on','Jconstant','on','Mass','M');
if keepall
t=linspace(0,T,N+1)';
[t,c]=solver('finode',t,c0,options,b,B,Phi);
c=c';
else
[t,c]=solver('finode',[0 T],c0,options,b,B,Phi);
c=c(end,:)';
end
V=Phi*c;
return
otherwise
error('Method option is invalid')
end
clear B M
% Invert the bases matrices to avoid repeated inversions
if tensor || american || barrier
for i=1:d
Phi{i}=inv(Phi{i});
% if the inverse matrix is not sparse, convert it to full
if nnz(Phi{i})/numel(Phi{i})>0.3,
Phi{i}=full(Phi{i});
end
end
end
% Initialize coefficient arrays
c=c0;
if keepall
C=zeros(n,N+1);
C(:,1)=c0;
end
% Iterate over time
for i=2:N+1
if hasdiv, c=a+A*c;
else c=A*c;
end
if mass
if tensor, c=ckronx(Phi,c);
else c=lusol(MLU,c);
end
end
if barrier || american
V=funeval(c,fspace,bases);
if barrier
V(Binda)=Ra;
V(Bindb)=Rb;
end
if american
if american<0
V=min(V0,V);
else
V=max(V0,V);
end
end
c=ckronx(Phi,V);
end
if keepall, C(:,i)=c; end
end
if keepall, c=C; end
V=funeval(c,fspace,bases);