rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
eop_example_script.m
Go to the documentation of this file.
1 function eop_example_script(step)
3 %
4 %% Step description:
5 % Attention: In some of the steps precomputed basis or other things are loaded!!
6 %
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
14 %
15 %% possible functionals l
16 %
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
21 %
22 % with G_0 being a subgrid of G the grid and P = (P^k)_{k=0}^M being the
23 % solution.
24 %
25 %% various offline-data included on the CD
26 %
27 % a model on the standard settings which contains both constants that are
28 % required for the various error estimators for basis generation and
29 % reduced siulation:
30 % - model_standard_settings_and_constants.mat
31 %
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
42 %
43 %
44 %% not included
45 %
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
49 % done in step 7.
50 %
51 
52 %% start of the script
53 % fields that can be changed by the user are commented with "choose" or
54 % "specify"
55 
56 text = sprintf('Starting step %d of eop_example_script at %s', step, datestr(now));
57 disp(text);
58 
59 % specify the grid and other elementary things of the model
60 params = [];
61 params.S1bar = 70;
62 params.S2bar = 70;
63 params.h1 = 0.5;
64 params.h2 = 0.5;
65 params.K = (params.S1bar + params.S2bar)/4;
66 model = european_option_pricing_model(params);
67 
68 switch step
69 
70  case 1 % primal or dual detailed simulation + plot. Can include detailed output in primal case.
71 
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.)
75 
76  if model.want_dual == 1 % primal case
77  model.compute_output_functional = 0; % choose whether to compute output (=1) or not (=0)
78  end
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';
83  end
84  % detailed solution
85  start = tic;
86  sim_data = detailed_simulation(model, model_data);
87  time_det_sim = toc(start);
88  disp(['detailed simulation took ' mat2str(time_det_sim) ' seconds'])
89  % plot
90  lin_evol_plot_sim_data(model,model_data,sim_data,params);
91 
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.
93 
94  clear
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.
96 
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)
100  end
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';
105  end
106 
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
110 
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)
114 
115  % generate reduced data, i.e. galerkinprojection on the reduced space
116  reduced_data = gen_reduced_data(model, detailed_data);
117 
118  model = set_mu(model,[0.05,0.4,2,0.5]); % choose the parameter
119  start = tic;
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
128  % plot output
129  figure;
130  plot(rb_sim_data.s)
131  title('output')
132  model.compute_output_functional = 1;
133  else
134  lin_evol_plot_sim_data(model,model_data,rb_sim_data,[]); % only plot simulation
135  end
136 
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.
140  if addition == 1
141  disp('addition choosen: detailed simulation gets computet...')
142  start = tic;
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
147  figure;
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)
150  axis tight
151  elseif strcmp(model.error_norm, 'energy') == 1 % error and estimator in space-time norm, i.e. just a number no sequence
152  error = eop_space_time_norm(rb_sim_data.U, sim_data.U, model_data.grid);
153  disp(['error = ' mat2str(error) ' estimator = ' mat2str(rb_sim_data.Delta)])
154  end
155  end
156 
157  case 3 % reduced basis generation
158 
159  mode = 1;
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).
162  if mode == 1
163  clear
164  load('model_standard_settings_and_constants.mat')
165  elseif mode == 2 % might be very expensive and RAM intense
166  model = eop_beta(model);
167  model = eop_C_L_I(model);
168  end
169 
170  model_data = gen_model_data(model); % high dimensional data (grid, etc.)
171 
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];
177  end
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)
182  end
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;
187  % RB generation
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
190 
191  case 4 % display a GUI to a loaded reduced basis. See eop_master_gui for more information on the GUI.
192 
193  clear;
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.
195  % prepare GUI input
196  sim_data = detailed_simulation(model, detailed_data);
197  reduced_data = gen_reduced_data(model, detailed_data);
198  % produce GUI
199  eop_master_gui(model, detailed_data, sim_data.U, reduced_data);
200 
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
202 
203  clear;
204  part = 1; % choose which part of the experiment you want to run (part2 can be very RAM intense)
205 
206  if part == 1 % part 1 - regards the EOC-computation and grid error (grid convergence)
207 
208  % specify starting model and number_of_refinementsteps
209  params = [];
210  params.number_of_refinementsteps = 6; % dont choose this too large
211  params.S1bar = 70;
212  params.S2bar = 70;
213  params.h1 = params.S1bar/2;
214  params.h2 = params.S2bar/2;
215  params.K = (params.S1bar + params.S2bar)/4;
216  model = european_option_pricing_model(params);
217  model_data = gen_model_data(model);
218  % compute grid error, see eop_gridconvergence for a detailed description
219  [grid_error, ~, EOC] = eop_gridconvergence(model, model_data, params);
220  % plot
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)]);
227 
228  elseif part == 2 % regarding the choice of the grid
229 
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
234  disp(i)
235  % sucsessivly build a larger model
236  params = [];
237  params.S1bar = i*10;
238  params.S2bar = i*10;
239  params.h1 = 0.5;
240  params.h2 = 0.5;
241  params.K = 20; % K is fix so that the various solutions are comparably
242  % generate and modify the model
243  model1 = european_option_pricing_model(params);
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
251  params.S1bar = i*15;
252  params.S2bar = i*15;
253  params.h1 = 0.5;
254  params.h2 = 0.5;
255  params.K = 20; % K is fix so that the various solutions are comparably
256  % generate and modify the model
257  model2 = european_option_pricing_model(params);
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)
267  for p = 1:model1.N2
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);
269  end
270  end
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);
273  end
274  semilogy(4:4+number_of_steps,rel_error,'r*')
275  axis tight
276  xlabel('i')
277  ylabel('relative error')
278 
279  end
280 
281 
282  case 6 % step showing that the solution is bounded uniformly in delta t - see Chapter 5.3 of the diploma thesis
283 
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
289  disp(i)
290  % refining \Delta t
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;
295  % computing ||P^M||
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);
298  end
299  % plot
300  semilogx(Delta_t,Norm,'r*')
301  axis tight
302  xlabel('\Delta t')
303  ylabel('|| P^M ||')
304 
305  case 7 % step computing the improved output and the improved output estimator. Also includes the unimproved output estimator.
306 
307  mode = 1;
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).
310  if mode == 1
311  clear
312  load('model_standard_settings_and_constants.mat')
313  elseif mode == 2 % might be very expensive and RAM intense
314  model = eop_beta(model);
315  model = eop_C_L_I(model);
316  end
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
336  detailed_data = lin_evol_primal_dual_gen_detailed_data(model, model_data); % dual and primal detailed_data
337  model.use_scm = 0; % doesn't work yet - but should be set here before gen_reduced_data (when SCM works)
338  % other fields regarding reduced_data generation are set in lin_evol_primal_dual_gen_reduced_data
339  reduced_data = lin_evol_primal_dual_gen_reduced_data(model, detailed_data); % dual and primal reduced_data
340  model = set_mu(model,[0.05,0.4,2,0.5]); % choose the parameter
341  rb_sim_data_improved = lin_evol_primal_dual_rb_simulation(model, reduced_data); % dual and primal rb_simulation
342  % some output
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
346  [~,l] = eop_fd_functionals(model, model_data);
347  if strcmp(model.name_output_functional, 'theta') == 1
348  non_improved_ouput_estimator = rb_sim_data_improved.Delta * 2*l;
349  else
350  non_improved_ouput_estimator = rb_sim_data_improved.Delta * l;
351  end
352  disp(['reduced improved output = ' mat2str(reduced_improved_output) ' improved outputestimator = ' mat2str(improved_output_estimator)...
353  ' unimproved outputestimator = ' mat2str(non_improved_ouput_estimator)])
354 
355  otherwise
356 
357  error('case unknown');
358 
359 end
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, ~)
Definition: eop_fd_norm.m:17
function [ model , beta_build ] = eop_beta(model)
[model, beta_build] = eop_beta(model)
Definition: eop_beta.m:17
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)
Definition: eop_C_L_I.m:17
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)