5 % Attention: In some of the steps precomputed basis or other things are loaded!!
7 % step 1: primal or dual detailed simulation with plot. Can include output in the primal
case.
8 % step 2: primal or dual reduced simulation with plot. No output here (is treated in step 7). Has an addition considering error and error estimator.
9 % step 3: generation and saving of primal (or dual) reduced basis.
10 % step 4: display of a GUI
for the reduced (and to some extend) the detailed simulation. Also contains error estimation.
11 % step 5: includes the experiments of Chapter 5.2 of the diploma thesis, i.e. stability analysis.
12 % step 6: includes the experiments of Chapter 5.3 of the diploma thesis, i.e. the boundedness of the solution
for Delta t -> 0.
13 % step 7: Treats the computation of the improved output and the improved output estimator
15 %% possible functionals l
17 % -
'average price': l(P^M;\mu) = 1/|\G_0|*\int_{\G_0}P^M(x)dx
18 % - 'theta': l(P^M;\mu) = -1/|\G_0|*int_{\G}\frac{\partial P^M}{\partial t}(x) dx
19 % -
'average dP/dS1': l(P^M;\mu) = 1/|\G_0|*\int_{\G_0}\frac{\parial P^M}{\partial S_1}(x)dx
20 % -
'average dP/dS2': l(P^M;\mu) = 1/|\G_0|*\int_{\G_0}\frac{\parial P^M}{\partial S_2}(x)dx
22 % with G_0 being a subgrid of G the grid and P = (P^k)_{k=0}^M being the
25 %% various offline-data included on the CD
27 % a model on the standard settings which contains both constants that are
28 % required
for the various error estimators
for basis generation and
30 % - model_standard_settings_and_constants.mat
32 % various files containing the 6 Reduced Basis (and corresponding model and
33 % model_data), described in Chapter 5.4. They use the error (err) or
34 % estimator (est) as error indicator and have different
35 % sizes (100, 200 or 300). The model always includes the empiric constants:
36 % - 70x70_h1h2_0.5_100_err.mat
37 % - 70x70_h1h2_0.5_100_est.mat
38 % - 70x70_h1h2_0.5_200_err.mat
39 % - 70x70_h1h2_0.5_200_est.mat
40 % - 70x70_h1h2_0.5_300_err.mat
41 % - 70x70_h1h2_0.5_300_est.mat
46 % a file containing a dual Reduced Basis (and corresponding model and
47 % model_data). But a dual Reduced Basis is only required in the case where
48 % the improved output and improved output estimator are desired and this is
52 %% start of the script
53 % fields that can be changed by the user are commented with "choose" or
59 % specify the grid and other elementary things of the model
65 params.K = (params.S1bar + params.S2bar)/4;
70 case 1 % primal or dual detailed simulation + plot. Can include detailed output in primal case.
72 model.want_dual = 0; % choose to compute a primal (model.want_dual = 0) or a dual (model.want_dual = 1) detailed simulation
73 model = set_mu(model,[0.05,0.4,0.5,2]); % choose the parameter
74 model_data = gen_model_data(model); % high dimensional data (grid, etc.)
76 if model.want_dual == 1 % primal case
77 model.compute_output_functional = 0; % choose whether to compute output (=1) or not (=0)
79 if model.compute_output_functional == 1 || model.want_dual == 1 % specify functional and subgrid G_0 in these cases
80 model.functional_subset_S1 = [15,15];
81 model.functional_subset_S2 = [15,15];
82 model.name_output_functional = 'average price';
86 sim_data = detailed_simulation(model, model_data);
87 time_det_sim = toc(start);
88 disp(['detailed simulation took ' mat2str(time_det_sim) ' seconds'])
90 lin_evol_plot_sim_data(model,model_data,sim_data,params);
92 case 2 % primal or dual reduced simulation + plot. Always load a precomputed reduced basis for the specific case. No improved output in this step.
95 load('70x70_h1h2_0.5_200_est.mat') % load a primal or dual reduced basis, computed with step 3 or an already existing one. Should contain the model, model_data and detailed_data.
97 model.want_dual = 0; % choose primal (=0) or dual (=1) case according to the data you just loaded
98 if model.want_dual == 0 % in primal case
99 model.compute_output_functional = 1; % choose whether to compute output (=1) or not (=0)
101 if model.compute_output_functional == 1 || model.want_dual == 1 % specify functional and subgrid in these cases
102 model.functional_subset_S1 = [15,15];
103 model.functional_subset_S2 = [15,15];
104 model.name_output_functional = 'average price';
107 % fields controlling the reduced_data generation - don't change them
108 model.use_scm = 0; % doesn't work yet.
109 model.want_improved_output = 0; % improved output is treated in step 7
111 model.error_norm = 'l2'; % choose 'l2' or 'energy' norm for the error estimator computed in the reduced simulation ('l2' is standard)
112 % for the energy error estimator a lower bound to the inf sup constant has to be available
113 % (is included in the included offline data via model.constant_LB)
115 % generate reduced data, i.e. galerkinprojection on the reduced space
116 reduced_data = gen_reduced_data(model, detailed_data);
118 model = set_mu(model,[0.05,0.4,2,0.5]); % choose the parameter
120 rb_sim_data = rb_simulation(model, reduced_data); % coefficients of reduced solution and error estimators
121 rb_sim_data = rb_reconstruction(model, detailed_data, rb_sim_data); % actual reduced solution
122 time_red_sim = toc(start);
123 disp(['reduced simulation took ' mat2str(time_red_sim) ' seconds'])
124 % neccessary since lin_evol_plot_sim_data can actually only plot detailed output
125 if model.compute_output_functional == 1 && model.want_dual == 0 % output only in primal case
126 model.compute_output_functional = 0;
127 lin_evol_plot_sim_data(model,model_data,rb_sim_data,[]); % plot the simulation
132 model.compute_output_functional = 1;
134 lin_evol_plot_sim_data(model,model_data,rb_sim_data,[]); % only plot simulation
137 % little addition - error computation and output of error and error
138 % estimator either as plot (if l2 norm) or simple numbers (if energy norm)
139 addition = 1; % addition is optional. Can be choosen as 1 or 0.
141 disp('addition choosen: detailed simulation gets computet...')
143 sim_data = detailed_simulation(model, model_data); % compute detailed simulation for error computation
144 time_det_sim = toc(start);
145 disp(['detailed simulation took ' mat2str(time_det_sim) ' seconds'])
146 if strcmp(model.error_norm, 'l2') == 1 % sequence of error and estimator
148 error =
eop_fd_norm(rb_sim_data.U, sim_data.U, model_data.grid);
149 semilogy(2:model.nt+1, error(2:end), 'r', 2:model.nt+1, rb_sim_data.Delta(2:end), 'b') % ignore the first timestep (estimator is always 0 there)
151 elseif strcmp(model.error_norm, 'energy') == 1 % error and estimator in space-time norm, i.e. just a number no sequence
153 disp(['error = ' mat2str(error) ' estimator = ' mat2str(rb_sim_data.Delta)])
157 case 3 % reduced basis generation
160 % choose to either load the standard model with the precomputed constants (mode = 1) or
161 % compute the constants for the error estimators for the current model (mode = 2).
164 load('model_standard_settings_and_constants.mat')
165 elseif mode == 2 % might be very expensive and RAM intense
170 model_data = gen_model_data(model); % high dimensional data (grid, etc.)
172 model.want_dual = 1; % choose to compute a dual (=1) or primal (=0) basis
173 if model.want_dual == 1 % for a dual basis specify the dual problem
174 model.name_output_functional = 'theta';
175 model.functional_subset_S1 = [15,15];
176 model.functional_subset_S2 = [15,15];
178 % specify some fields controlling the basis generation.
179 model.RB_error_indicator = 'estimator';
180 if strcmp(model.RB_error_indicator, 'estimator') == 1
181 model.error_norm = 'energy'; % 'l2' or 'energy' norm for the error estimator computed in the reduced simulation ('l2' is standard)
183 model.RB_numintervals = 5*ones(length(model.mu_ranges),1); % fineness of the training set
184 model.RB_stop_epsilon = eps;
185 model.RB_stop_Nmax = 10;
186 model.RB_stop_timeout = 120*3600;
188 detailed_data = gen_detailed_data(model, model_data);
189 save('how_you_want_to_name_your_basis.mat','model','model_data','detailed_data') % the
string specifys the filename
191 case 4 % display a GUI to a loaded reduced basis. See eop_master_gui for more information on the GUI.
194 load('70x70_h1h2_0.5_200_est.mat') % load a primal or dual reduced basis, computed with step 3 or an already existing one. Should contain the model, model_data and detailed_data.
196 sim_data = detailed_simulation(model, detailed_data);
197 reduced_data = gen_reduced_data(model, detailed_data);
199 eop_master_gui(model, detailed_data, sim_data.U, reduced_data);
201 case 5 % this step shows the grid-convergence, EOC (exponential order of convergence) and choice of grid - see Chapter 5.2 of the diploma thesis
204 part = 1; % choose which part of the experiment you want to run (part2 can be very RAM intense)
206 if part == 1 % part 1 - regards the EOC-computation and grid error (grid convergence)
208 % specify starting model and number_of_refinementsteps
210 params.number_of_refinementsteps = 6; % dont choose this too large
213 params.h1 = params.S1bar/2;
214 params.h2 = params.S2bar/2;
215 params.K = (params.S1bar + params.S2bar)/4;
217 model_data = gen_model_data(model);
221 semilogy(1:params.number_of_refinementsteps, grid_error, 'r');
222 set(gcf,'Color','white');
223 xlabel('number of refinement steps');
224 set(gca,'XTick',1:params.number_of_refinementsteps);
225 ylabel('grid error');
226 disp(['EOC-sequence:'
char(10) mat2str(EOC)]);
228 elseif part == 2 % regarding the choice of the grid
230 number_of_steps = 6; % choose number of procedure steps - dont choose this too large
231 mu = [0.05,0.4,0.5,0.5]; % choose the parameter for the experiment
232 rel_error = zeros(number_of_steps,1);
233 for i = 4:4+number_of_steps
235 % sucsessivly build a larger model
241 params.K = 20; % K is fix so that the various solutions are comparably
242 % generate and modify the model
244 model1.save_time_indices = model1.nt; % only save the solution for the last point in time
245 model1 = set_mu(model1,mu);
246 model1.init_value_function = 'standard';
247 model_data1 = gen_model_data(model1);
248 % detailed simulation on the coarse grid
249 sim_data_coarse = detailed_simulation(model1, model_data1);
250 % now build the model on the refined grid
255 params.K = 20; % K is fix so that the various solutions are comparably
256 % generate and modify the model
258 model2.save_time_indices = model2.nt; % only save the solution for the last point in time
259 model2 = set_mu(model2,mu);
260 model2.init_value_function = 'standard';
261 model_data2 = gen_model_data(model2);
262 % detailed simulation on the refined grid
263 sim_data_refined = detailed_simulation(model2, model_data2);
264 % cutoff of the refined data
265 sim_data_refined_cutoff = zeros(size(sim_data_coarse.U));
266 for o = 1:size(sim_data_refined_cutoff,2)
268 sim_data_refined_cutoff((p-1)*model1.N1+1 : p*model1.N1,o) = sim_data_refined.U((p-1)*model2.N1+1 : (p-1)*model2.N1 + model1.N1,o);
271 % compute relative error
272 rel_error(i-3) =
eop_fd_norm(sim_data_coarse.U, sim_data_refined_cutoff, model_data1.grid)/
eop_fd_norm(sim_data_coarse.U, zeros(size(sim_data_coarse.U)), model_data1.grid);
274 semilogy(4:4+number_of_steps,rel_error,'r*')
277 ylabel('relative error')
282 case 6 % step showing that the solution is bounded uniformly in delta t - see Chapter 5.3 of the diploma thesis
284 number_of_steps = 10; % choose the number of steps, meaning i = 1:number_of_steps
285 model_data = gen_model_data(model);
286 Delta_t = zeros(number_of_steps,1);
287 Norm = zeros(number_of_steps,1);
288 for i = 1:number_of_steps
291 model.deltaT = 2^(-i);
292 Delta_t(i) = model.deltaT;
293 model.nt = model.T/model.deltaT;
294 model.save_time_indices = model.nt;
296 sim_data = detailed_simulation(model, model_data);
297 Norm(i) =
eop_fd_norm(sim_data.U(:,end), zeros(size(sim_data.U(:,end))), model_data.grid);
300 semilogx(Delta_t,Norm,'r*')
305 case 7 % step computing the improved output and the improved output estimator. Also includes the unimproved output estimator.
308 % choose to either load the standard model with the precomputed constants (mode = 1) or
309 % compute the constants for the error estimators for the current model (mode = 2).
312 load('model_standard_settings_and_constants.mat')
313 elseif mode == 2 % might be very expensive and RAM intense
317 model_data = gen_model_data(model); % high dimensional data (grid, etc.)
318 % specify the dual problem
319 model.name_output_functional = 'theta';
320 model.functional_subset_S1 = [15,15];
321 model.functional_subset_S2 = [15,15];
322 % RB settings (primal and dual) - can be changed
323 model.RB_stop_epsilon = eps;
324 model.RB_stop_Nmax = 200;
325 model.RB_stop_timeout = 3600*3600;
326 model.dual_RB_stop_epsilon = eps;
327 model.dual_RB_stop_Nmax = 200;
328 model.dual_RB_stop_timeout = 3600*3600;
329 % dont change the following 3 fields unless you know what you are doing
330 model.RB_numintervals = 5*ones(length(model.mu_ranges),1);
331 model.RB_error_indicator = 'estimator';
332 model.error_norm = 'energy'; % use energy estimator for basis generation since energy error estimator is used later for the improved output estimator
333 % compute primal/dual detailed/reduced data so that the improved output and output estimator can be computed via
335 % See the respective m-files for more information
337 model.use_scm = 0; % doesn't work yet - but should be set here before gen_reduced_data (when SCM works)
340 model = set_mu(model,[0.05,0.4,2,0.5]); % choose the parameter
343 reduced_improved_output = rb_sim_data_improved.s;
344 improved_output_estimator = rb_sim_data_improved.Delta_s;
345 % need the continuousity constant of the functional for the non_improved_ouput_estimator
347 if strcmp(model.name_output_functional, 'theta') == 1
348 non_improved_ouput_estimator = rb_sim_data_improved.Delta * 2*l;
350 non_improved_ouput_estimator = rb_sim_data_improved.Delta * l;
352 disp(['reduced improved output = ' mat2str(reduced_improved_output) ' improved outputestimator = ' mat2str(improved_output_estimator)...
353 ' unimproved outputestimator = ' mat2str(non_improved_ouput_estimator)])
357 error('case unknown');
function reduced_data = lin_evol_primal_dual_gen_reduced_data(model, detailed_data)
reduced_data = lin_evol_primal_dual_gen_reduced_data(model, detailed_data)
function [ v , l ] = eop_fd_functionals(model, model_data)
[v, l] = eop_fd_functionals(model, model_data)
function model = european_option_pricing_model(params)
model = european_option_pricing_model(params)
function norm = eop_fd_norm(fd_function1, fd_function2,gridbase grid, unused1)
norm = eop_fd_norm(fd_function1, fd_function2, grid, ~)
function [ model , beta_build ] = eop_beta(model)
[model, beta_build] = eop_beta(model)
function simulation_data = lin_evol_primal_dual_rb_simulation(model, reduced_data)
simulation_data = lin_evol_primal_dual_rb_simulation(model, reduced_data)
function model = eop_C_L_I(model)
model = eop_C_L_I(model)
function norm = eop_space_time_norm(fd_function1, fd_function2,gridbase grid, unused1)
norm = eop_space_time_norm(fd_function1, fd_function2, grid, ~)
function eop_example_script(step)
eop_example_script(step)
function [ grid_error , h_sequence , EOC ] = eop_gridconvergence(model, model_data, params)
[error, h_sequence, EOC] = eop_gridconvergence(model, model_data, params)
function detailed_data = lin_evol_primal_dual_gen_detailed_data(model, model_data)
detailed_data = lin_evol_primal_dual_gen_detailed_data(model, model_data)