eop_beta.m

Go to the documentation of this file.

function [model, beta_build] = eop_beta(model)
%[model, beta_build] = eop_beta(model)
%
% function computing the inf-sup constant of the PG-BLF-Matrix of the
% european_option_pricing_model for each mu in a training set. Output is
% the minimum of those inf-sup constants - the empirical inf-sup constant.
% Uses the huge sparse matrix and therefor needs a lot of RAM. Could also
% be done with AFUN for eigs but this approach was sufficient for one run
% to get an empirical constant.

% Dominik Garmatter 17.07.2012

model_data = gen_model_data(model);
K = model.get_inner_product_matrix(model_data);
% get the components
model.decomp_mode = 1;
[L_I_comp, L_E_comp] = model.operators_ptr(model, model_data);
% M*H is the size of the PG-BLF-Matrix
H = size(L_I_comp{1},1);
M = model.nt+1;
% build Mtrain as a first parameter and many more random parameters
build1 = zeros(size(model.mu_ranges,2),1);
build2 = zeros(size(model.mu_ranges,2),1);
for i = 1:size(model.mu_ranges,2)
    build1(i,1) = model.mu_ranges{i}(1);
    build2(i,1) = model.mu_ranges{i}(2);
end
Mtrain1 = [0.05,0.4,2,0.5]';
Mtrain2 = unifrnd(repmat(build1,1,50),repmat(build2,1,50));
Mtrain = [Mtrain1, Mtrain2];

beta_build = zeros(size(Mtrain,2),1);
beta_build2 = zeros(size(Mtrain,2),1);

kron_1 = sparse(diag([1;zeros(model.nt,1)]));
kron_2 = sparse(diag([0;ones(model.nt,1)]));
kron_3 = sparse(diag(-ones(model.nt,1),-1));

for k = 1:size(Mtrain,2)
    disp(k);
    model = model.set_mu(model,Mtrain(:,k));
    % compute the current coefficients
    old_decomp_mode = model.decomp_mode;
    model.decomp_mode = 2;
    [L_I_coeffs, L_E_coeffs] = model.operators_ptr(model);
    model.decomp_mode = old_decomp_mode;
    % and assemble the matrices
    L_I = lincomb_sequence(L_I_comp, L_I_coeffs);
    L_E = lincomb_sequence(L_E_comp, L_E_coeffs);
    % compute the smallest eigenvalue of B'*K*B and take the sqareroot of it
    Matrix = kron(kron_1, speye(H)) + kron(kron_2, L_I) + kron(kron_3, L_E);
    try %first try to solve the EWP with default tolerance eps
      beta_build(k) = sqrt(eigs(Matrix'*kron(speye(M),K)*Matrix,1,'SM'));
    catch err %if no eigenvalues to sufficient accuracy was found reduce the tolerance and try it again
      if (strcmp(err.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) 
          opts.tol = 1e-14;
          no_solution = 1;
          while no_solution
            disp(['I reduced the tolerance of an eigenvalueproblem to ' mat2str(opts.tol)]);
            try
              beta_build2(k) = sqrt(eigs(Matrix'*kron(speye(M),K)*Matrix,1,'SM',opts));
              no_solution  = 0; %end when there is no error
            catch err2 %when it still fails reduce the tolerance and try again!
              no_solution = 1;
              if (strcmp(err2.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14'))  
                opts.tol = 100 * opts.tol;
              end
              if opts.tol > 1e-3 %if the tolerance is reduced until a certain value, give an error since the tolerance is too soft!
                  error(['I reduced the tolerance to ' mat2str(opts.tol) ' but I still could not find a solution. Reconsider your Problem']);
              end
            end
          end
      else %if there is another error rethrow
          rethrow(err)
      end
    end
end

model.beta_LB = min(beta_build); % take the minimum