european_option_pricing_model.m

Go to the documentation of this file.

function model = european_option_pricing_model(params)
%model = european_option_pricing_model(params)
%
% The model of my diploma thesis for the solution of the two-dimensional
% Black-Scholes-Equation.
%
% See Chapter 2.3 and 4.1 for a precise description of the problem.
%
% Dominik Garmatter 20.08 2012

%%                     GENERAL SETTINGS

model.rb_problem_type           = 'lin_evol';
model.verbose                   = 5; 
model.affinely_decomposed       = true;
model.compute_output_functional = 0;
model.debug                     = 0;

%%                      MU RELATED

model.mu_names  = {'r','rho','sigma1','sigma2'}; % r = the interest rate, roh the correlation coefficient, sigma 1 und 2 the volatilities
model.mu_ranges = {[1e-10 0.3],[-0.9 0.9],[1e-10 2],[1e-10 2]}; % the parameter space
% some initial values
model.r         = 0.05;
model.rho       = 0.4;
model.sigma1    = 0.5; 
model.sigma2    = 0.5;
% functions to set or get a mu
model.set_mu    = @set_mu_default;
model.get_mu    = @get_mu_default;

%%                      TIME RELATED

model.T                  = 1;                % date of maturity
model.nt                 = 20;               % number of timeintervals
model.deltaT             = model.T/model.nt;
model.data_const_in_time = 1;
model.set_time           = @set_time_default;

%%                  GRID RELATED SETTINGS

model.gridtype = 'rectgrid';

if nargin == 0  
    params = [];
end

% grid construction by using a params-structure, which can be given as an
% input to the model. The structure contains the following fields:
%- S1bar -> size of the grid in S1-direction (x-direction)
%- S2bar -> size of the grid in S2-direction (y-direction)
%- h1    -> step-width in S1-direction
%- h2    -> step-width in S2-direction
%- K     -> the Strike (not grid relevant and normally choosen as (params.S1bar + params.S2bar)/4)

if ~isfield(params,'S1bar')
    params.S1bar     = 70;
end
params.xrange        = [0,params.S1bar];

if ~isfield(params,'S2bar')
    params.S2bar     = 70;
end
params.yrange        = [0,params.S2bar];

if ~isfield(params,'h1')
    params.h1        = 0.5;
end
params.xnumintervals = params.S1bar/params.h1;

if ~isfield(params,'h2')
    params.h2        = 0.5;
end
params.ynumintervals = params.S2bar/params.h2;

if ~isfield(params,'K')
    params.K         = (params.S1bar + params.S2bar)/4;
end

model                = structcpy(model,params);
model.N1             = model.xnumintervals + 1; % number of gridpoints in S1-direction
model.N2             = model.ynumintervals + 1; % number of gridpoints in S2-direction

%%                         RB RELATED

model.RB_generation_mode          = 'greedy_uniform_fixed';
model.RB_extension_algorithm      = @RB_extension_PCA_fixspace;
model.set_rb_in_detailed_data     = @(detailed_data,RB) ...
                                    setfield(detailed_data, 'RB', RB);
model.get_rb_from_detailed_data   = @(detailed_data) detailed_data.RB;
model.get_rb_size                 = @get_rb_size_default;
model.get_dofs_from_sim_data      = @(sim_data) sim_data.U;
model.rb_init_values              = @rb_init_values_primal_dual; % includes primal/dual difference
model.rb_init_data_basis          = @RB_init_data_basis;
model.RB_numintervals             = 5*ones(length(model.mu_ranges),1); % fineness of Mtrain

model.RB_error_indicator          = 'estimator';
% fields required in case of model.RB_error_indicator = 'error'
model.error_algorithm             = @eop_fd_norm; % discrete L2-norm
model.RB_detailed_train_savepath  = 'detailed_sim_from_basisgen';
model.save_detailed_simulations   = @save_detailed_simulations;
% fields required in case of model.RB_error_indicator = 'estimator'
model.error_norm                  = 'l2'; % optional 'energy'
model.get_estimator_from_sim_data = @(sim_data) sim_data.Delta(end);  

% RB stopping conditions                                        
model.RB_stop_epsilon             = eps;
model.RB_stop_Nmax                = 100;
model.RB_stop_timeout             = 3600*3600;
% stopping conditions for the dual RB
model.dual_RB_stop_epsilon        = eps;
model.dual_RB_stop_Nmax           = 100;
model.dual_RB_stop_timeout        = 3600*3600;

% detailed_data get saved every time - this specifies the filename they get saved in
model.detailedfname               = 'eop_detailed_data';   

% following fields get removed temporary so filecaching during the basisgeneration works
model.filecache_ignore_fields_in_model = {'N','Nmax'};
model.filecache_ignore_fields_in_detailed_data = {'RB_info'};


%%                         VARIOUS

model.gen_model_data              = @lin_evol_gen_model_data; % generating the mass matrix and the grid
model.mass_matrix                 = @eop_mass_matrix;
model.inner_product               = @(model, model_data, coeffvector1, coeffvector2)...  % this is the discrete L2-inner-product:
                                    coeffvector1'*model_data.W*coeffvector2;             %<u,v> = 1/H \sum_{x\in G}{u(x)*v(x)}, with G being the grid and H = #G
model.get_inner_product_matrix    = @(detailed_data) detailed_data.W;
model.detailed_simulation         = @lin_evol_detailed_simulation_primal_dual; % includes primal/dual difference     
model.gen_detailed_data           = @lin_evol_gen_detailed_data;
model.orthonormalize              = @model_orthonormalize_gram_schmidt;
model.gen_reduced_data            = @lin_evol_gen_reduced_data_primal_dual; % includes primal/dual difference  
model.rb_operators                = @lin_evol_rb_operators_primal_dual; % includes primal/dual difference
model.rb_simulation               = @lin_evol_rb_simulation_primal_dual; % includes primal/dual difference
model.rb_reconstruction           = @rb_reconstruction_default; 
model.reduced_data_subset         = @lin_evol_reduced_data_subset;      
  
%%                         FD & Functionals  

model.init_values_algorithm    = @eop_init_values; % includes 2 different functions
model.init_value_function      = 'standard'; % see "eop_init_values"
model.operators_ptr            = @eop_fd_operators; % these compute L_I,L_E and b
model.operators_output         = @eop_fd_functionals; % includes various output functionals
model.name_output_functional   = 'average price'; % see "eop_fd_functionals"

% specify G_0 the subgrid on which the functionals are evaluated
model.functional_subset_S1 = model.xrange; 
model.functional_subset_S2 = model.yrange;

%%                         PLOT RELATED

model.plot                  = @plot_vertex_data;
model.plot_sim_data         = @lin_evol_plot_sim_data;

%%                         NORM BOUNDS (for error estimators)

% for the empiric constants use eop_C_L_I and eop_beta. For the SCM
% (model.use_scm = 1) see scm_offline and lin_evol_gen_reduced_data_primal_dual.
% But right now SCM is not working -> use the empiric constants.
model.L_I_inv_norm_bound = 1.000431891324076; % empiric constant for standard settings. Use eop_C_L_I to compute an empiric constant for other settings.
model.L_E_norm_bound = 1; % L_E = Id
model.constant_LB = 6.448207658122034e-005; % empiric LowerBound for standard settings. Use eop_beta to compute an empiric LowerBound for other settings.

end


function W = eop_mass_matrix(~,grid,~)
% function generating the mass matrix to the discrete Basis 
% \psi_i(x)= 1 , if x=x_i , 0 else
% and the discrete inner product <u,v> = 1/H \sum_{x\in G}{u(x)*v(x)}, 
% with G being the grid and H = #G.

W = speye(grid.nvertices)*1/grid.nvertices;
end