% calibrate_volatility_break_2.m
% Aug 3, 2014, MWW
% Calibrate volatility experiment 2 .. break in variance of I(1) component
% ----------- 
clear all;
nmc = 100000;

T = 260; 
scl = 0.707;   % Chosen so that initial theta = 0.5;



p = 1/260;
delta_vec = [0.814 1.612]';

p = 1/40;
delta_vec = [0.407 0.806]';

for i = 1:size(delta_vec,1);
 delta = delta_vec(i);
 s = rand(T,nmc) < p;
 tmp = rand(T,nmc);
 sgn_ind = (tmp < 0.5) - (tmp >= 0.5);
 del_mu = (s.*sgn_ind)*delta;
 mu_T = sum(del_mu)';
 sig_T = exp(mu_T);
 sig_T_scl = scl*sig_T;

 pct = [0.25;0.50;0.75];
 rslt_pctile = pctile(sig_T,pct);
 rslt_iqr = rslt_pctile(3)-rslt_pctile(1);
 fprintf('delta = %7.4f \n',delta);
 fprintf('75th percentile = %6.3f \n',rslt_pctile(3));
 fprintf('50th percentile = %6.3f \n',rslt_pctile(2));
 fprintf('25th percentile = %6.3f \n',rslt_pctile(1));
 fprintf('IQR = %6.3f \n',rslt_iqr);

 % Convert to theta values 
 var_eps = sig_T_scl.^2;
 var_eta = ones(nmc,1);
 % MA Coefficient and Standard Deviation of MA innovation 
 lam0=var_eps+2*var_eta;
 rho1=-var_eta./lam0;
 theta=-(ones(size(rho1))-sqrt(ones(size(rho1))-4*rho1.^2))./(2*rho1);
 var_a=lam0./(ones(size(rho1))+theta.^2);
 sd_a=sqrt(var_a);
 pct = [0.25;0.50;0.75];
 rslt_pctile = pctile(theta,pct);
 rslt_iqr = rslt_pctile(3)-rslt_pctile(1);
 fprintf('75th percentile (theta) = %6.3f \n',rslt_pctile(3));
 fprintf('50th percentile (theta) = %6.3f \n',rslt_pctile(2));
 fprintf('25th percentile (theta) = %6.3f \n',rslt_pctile(1));
 fprintf('IQR  (theta) = %6.3f \n',rslt_iqr);
 fprintf('\n');
end;
 
