Home > code > main > wasc_monte_carlo.m

wasc_monte_carlo

PURPOSE ^

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 cd('..');
0002 boot;
0003 cd(cgmm_config.directories.main);
0004 
0005 % Init parameters
0006 n = 2
0007 dt = 1/250;
0008 
0009 % load WASC parameters estimated from the corresponding time series
0010 % to be used as 'true' parameters for simulation
0011 load(cgmm_config.estimates.wasc);
0012 
0013 % rename estimates to make clear that they play the role of true parameters here
0014 mu = mu_cgmm
0015 Sigma_0 = Sigma_0_cgmm
0016 M = M_cgmm
0017 Q = Q_cgmm
0018 rho = rho_cgmm
0019 beta = round(beta_cgmm) % integer beta required for simulation
0020 
0021 % prepare simulation parameters
0022 S_0 = [100 100];
0023 y_0 = log(S_0);
0024 time_steps = cgmm_config.monte_carlo.time_steps/dt;
0025 
0026 % encode true parameters into flat parameter vector
0027 [theta_flat_0, decode] = encode_wasc_param(mu, Sigma_0, M, Q, rho, beta);
0028 % prepare characteristic function call with the encoded parameters
0029 cf = @(omega, th, y_t, tau) cf_wasc_v_theta(decode, omega, th, y_t, tau);
0030 % prepare options for optimization routine
0031 options = optimset('Display', 'iter' ...
0032                   , 'Algorithm', 'interior-point');
0033 % contraint for all parameters of +-10% around heuristic estimate
0034 lb = theta_flat_0 - abs(theta_flat_0)*0.1;
0035 ub = theta_flat_0 + abs(theta_flat_0)*0.1;
0036 %lb(6) = -max(abs(theta_flat_0(5:8)))*0.1; % lower left Q; heurist estimate = 0
0037 %ub(6) = max(abs(theta_flat_0(5:8)))*0.1; % lower left Q; heurist estimate = 0
0038 % except for correlation (where absolute constraints are more feasible)
0039 %lb(end-2:end-1) = max(theta_flat_0(end-2:end-1)-0.6, -1);
0040 %ub(end-2:end-1) = min(theta_flat_0(end-2:end-1)+0.6, 1);
0041 
0042 simulation_runs = cgmm_config.monte_carlo.simulation_runs
0043 if exist(cgmm_config.monte_carlo.wasc)
0044   load(cgmm_config.monte_carlo.wasc);
0045   % loaded file contains cell array 'estimates'
0046   first_run = size(estimates{1},1) + 1
0047 else
0048   estimates = {}; % store a cell of estimates for each time step
0049   % each estimates{k} is a simulation_runs x number of params matrix
0050   for k = 1:length(time_steps)
0051     estimates{k} = theta_flat_0;
0052   end
0053   first_run = 1
0054 end
0055 
0056 for run = first_run:simulation_runs
0057   disp(strcat('Simulation run: ',num2str(run),'/',num2str(simulation_runs)));
0058   for k = 1:length(time_steps)
0059     t = 0:dt:(time_steps(k)*dt);
0060     % simulate time series
0061     y = sim_wasc_2d(y_0, mu, Sigma_0, M, Q, rho, beta, t);
0062     % perform parameter estimation
0063     tic;
0064     [theta_flat_cgmm, theta_flat_first] = cgmm(y, dt, cf, theta_flat_0 ...
0065                                             , cgmm_config.cgmm.grid_min ...
0066                                             , cgmm_config.cgmm.grid_max ...
0067                                             , cgmm_config.cgmm.grid_res ...
0068                                             , lb, ub, options);
0069     %theta_flat_cgmm = theta_flat_first = theta_flat_0;
0070     toc;
0071     estimates{k} = [estimates{k}; theta_flat_first];
0072   end
0073   % save monte carlo estimates in each iteration
0074   save(cgmm_config.monte_carlo.wasc, 'estimates', 'time_steps');
0075 end
0076

Generated on Mon 29-Apr-2013 19:29:13 by m2html © 2005