%% CHJ GARCH model simulations (Table IV, part 1)
%
% Simulates the Christoffersen-Heston-Jacobs (2013) GARCH model to compute
% phi and gamma bounds from option-implied beliefs. True phi is estimated
% by conditioning on initial volatility quantiles
%
% Input:  None (self-contained simulation)
% Output: Output/table4_part1.txt, ../FinalTables/table4_part1.txt
%         Output/Simulations/garch_output.mat

clear;
close all;
addpath(genpath(pwd));
resimulate = 1; % set to 1 to rerun sims
resimulate_short = 1; % set to 1 to just run 200 simulations quickly; set to 0 to run the full set of 10,000 sims
estimate_true = 1; % set to 1 to redraw S and h to estimate true probabilities (keep equal to 1)
warning('off','MATLAB:integral:NonFiniteValue');

% Simulation parameters
n = 10000*(1-resimulate_short) + 200*resimulate_short; % number of sims
n_full = n;
n_tot = 4; % split into four sub-simulation exercises to allow save (increment n_i each time)
dpy = 252; % trading days per year
burn_yr = 10; % burn-in period in years
burn = burn_yr*dpy; % convert to daily
sample_yr = 23*10; % number of years in true sample (true sample is 23 years, but want more since will do non-overlapping contracts)
sample = sample_yr*dpy;
sim_T = burn+sample; % simulation length
seednum = 3027; % arbitrary seed
rng(seednum,"twister"); % rng seed
T = 21*3; % 3-month contracts (21 days/month)
exp_dates = zeros(sim_T + 1,1);
exp_dates(burn+1:T:end) = 1; % vector selecting expiration dates
exp_count = sum(exp_dates>0);
exp_inds = find(exp_dates>0);

% Parameters from CHJ (2013)
omega = 0;
alpha = 8.887e-7;
beta = 0.756;
gamma = 515.57;
mu = 1.594; %  model-implied risk aversion (equity premium = mu*sigma^2)
chi_mult = 1.2638; % (1-2*alpha*chi)^-1
o_star = omega*chi_mult;
a_star = alpha*chi_mult^2;
b_star = beta;
phi = -(mu-1/2+gamma)/chi_mult + gamma - 1/2;
g_star = gamma-phi;
h_0 = (omega+alpha)/(1-beta-alpha*gamma^2); % initialize at long-run mean
h_0_star = (o_star+a_star)/(1-beta-a_star*g_star^2); % same
S_0 = 100; % initialize at 100
rf = 0.05/dpy;
param = [a_star b_star g_star o_star];

%% Simulations

if resimulate==1

    % Setup for spot draws
    S_sims = zeros(sim_T + 1,n); % each col. has one sim (burn-in plus true plus init)
    S_sims(1,:) = S_0;
    h_sims = zeros(sim_T + 1,n);
    h_sims(1,:) = h_0;
    z_sims = [zeros(1,n);randn(sim_T,n)];
    
    % RN volatility and price paths
    for t=2:sim_T+1
        h_sims(t,:) = omega + beta*h_sims(t-1,:) + alpha*(z_sims(t-1,:)-gamma*sqrt(h_sims(t-1,:))).^2;
        S_sims(t,:) = exp(log(S_sims(t-1,:)) + rf + (mu-1/2)*h_sims(t,:) + sqrt(h_sims(t,:)).*z_sims(t,:));
    end

    S_sims_full = S_sims;
    h_sims_full = h_sims;

    for n_i = 1:4
        disp(n_i);
        
        % Cut up into sub-simulation (do this here to avoid having to reseed)
        n_inds = (n_i-1)*n_full/n_tot+1:n_i*n_full/n_tot;
        S_sims = S_sims_full(:,n_inds);
        h_sims = h_sims_full(:,n_inds);
        n = n_full/n_tot;
        
        % Diagnostic plots
        diag_plots = 0; % expand for loop below for plots code
        for dd=1
            if diag_plots==1
                figure(1)
                plot(S_sims(:,10));
                figure(2)
                plot(sqrt(h_sims(:,10)));
                figure(3)
                plot(mean(S_sims(burn+2:10000,:),2));
                figure(4)
                plot(mean(sqrt(h_sims(burn+2:10000,:)),2));
                
                S_use = S_sims(5000,10);
                K = S_use*(0.75:0.01:1.25)';
                S = S_use;
                h_use = h_sims(5001,10);
                tau = 60;
                
                [rn,c] = RN_CHJ(S,K,h_use,rf,tau,param,1);
                figure(5)
                plot(K,c);
                figure(6)
                plot(K,rn);
                figure(7)
                plot(K(2:end),diff(rn));
            end
        end
        
        % Setup for option price calcs
        realized_ret = log(S_sims(exp_inds(2:end),:)./S_sims(exp_inds(1:end-1),:))-rf*T; % realized excess returns over option contracts
        max_return = 1; % min/max excess return capped at 100% for numerical integration
        bin_edges = [-max_return -0.2:0.05:0.2 max_return]; % discretize into 5-ppt ranges
        realized_bins = discretize(realized_ret,bin_edges);
        realized_ret_trunc = realized_ret; % truncate just for calculating strikes at which to calculate RN distribution
        realized_ret_trunc(realized_ret_trunc<=-0.25) = -0.225;
        realized_ret_trunc(realized_ret_trunc>=0.25) = 0.225;
        
        K_min_pct = discretize(realized_ret_trunc-0.05,bin_edges,bin_edges(1:end-1)); % will calculate at four strikes bracketing the three bins centered at the realized one
        K_2_pct = discretize(realized_ret_trunc,bin_edges,bin_edges(1:end-1));
        K_3_pct = discretize(realized_ret_trunc,bin_edges,bin_edges(2:end));
        K_max_pct = discretize(realized_ret_trunc+0.05,bin_edges,bin_edges(2:end));
        K_min_pct(K_2_pct==-max_return) = NaN; % only 3 strikes needed if return is <-20% or >20%
        K_max_pct(K_3_pct==max_return) = NaN;
        start_level = S_sims(exp_inds(1:end-1),:);
        start_array = repmat(start_level,[1 1 4]);
        strikes_pre = start_array.*exp(rf*T+cat(3,K_min_pct,K_2_pct,K_3_pct,K_max_pct)); % convert strikes from percent excess returns to terminal prices given starting prices
        strikes = repmat(strikes_pre,[1 1 1 T]);
        rn_cdf = NaN(exp_count-1,n,4,T); % dims: 1 for contract #, 2 for simulation #, 3 for strike, 4 for trading day within contract
        dtm = repmat(shiftdim(T:-1:1,-2),[exp_count-1,n,4,1]);
        rn_bin_num = cat(3,realized_bins-1,realized_bins,realized_bins+1); % keep track of bin numbers
        rn_bin_num(rn_bin_num==0|rn_bin_num>max(max(realized_bins))) = NaN;

        % Run through by dtm to calculate RN CDFs
        for j=T:-1:1
            disp(j);
            dtm_j = T-j+1;
            dt_inds = exp_inds(2:end)-dtm_j;
            prices_j = S_sims(dt_inds,:);
            price_slice = repmat(prices_j,[1 1 4]);
            vols_j = h_sims(dt_inds,:);
            vol_slice = repmat(vols_j,[1 1 4]);
            dtm_slice = dtm(:,:,:,j);
            strike_slice = strikes(:,:,:,j);
            [rn_out,~] = RN_CHJ(price_slice(:),strike_slice(:),vol_slice(:),rf,dtm_slice(:),param,0);
            rn_cdf(:,:,:,j) = reshape(rn_out,[exp_count-1,n,4,1]);
        end
        
        rn_probs = diff(rn_cdf,1,3); % 3rd dim is bin (low/med/high, where med is realized)
        realization = repmat(shiftdim([0 1 0],-1),[exp_count-1 n 1 1]); % for last day
        rn_probs = cat(4,rn_probs,realization);
        rn_probs(rn_probs<0)=0;
        rn_probs(rn_probs>1)=1;
        
        rn_binary = rn_probs(:,:,1:2,:)./(rn_probs(:,:,1:2,:)+rn_probs(:,:,2:3,:));
        m_daily = diff(rn_binary,1,4).^2;
        u_daily = rn_binary.*(1-rn_binary);
        r_daily = -diff(u_daily,1,4);
        X_daily = m_daily-r_daily;
        X_contract = sum(X_daily,4,"omitnan");
        u_0 = u_daily(:,:,:,1);
        pi_0 = rn_binary(:,:,:,1);
        
        file_init = "Output/Simulations/garch_sims_%d.mat";
        file_nm = sprintf(file_init,n_i);
        save(file_nm,"X_contract","pi_0","u_0","rn_bin_num","-v7.3");

    end

    % Reset to correct number of sims
    n = n_full;
end

%% Load and combine sub-simulations
X_contract_full = NaN(exp_count-1,n,2);
pi_0_full = NaN(exp_count-1,n,2);
u_0_full = NaN(exp_count-1,n,2);
rn_bin_num_full = NaN(exp_count-1,n,3);
file_init = "Output/Simulations/garch_sims_%d.mat";
for j=1:4
    file_nm = sprintf(file_init,j);
    sims_imp = load(file_nm);
    n_inds = (j-1)*n/n_tot+1:j*n/n_tot;
    X_contract_full(:,n_inds,:) = sims_imp.X_contract;
    pi_0_full(:,n_inds,:) = sims_imp.pi_0;
    u_0_full(:,n_inds,:) = sims_imp.u_0;
    rn_bin_num_full(:,n_inds,:) = sims_imp.rn_bin_num;    
end

rn_bin_pair = rn_bin_num_full(:,:,1:2);
rn_bin_pair(isnan(u_0_full)) = NaN;
rn_bin_use = ~isnan(rn_bin_pair) & rn_bin_pair>=2 & rn_bin_pair<=8;

X_pct = mean(X_contract_full(X_contract_full~=0))/mean(u_0_full(~isnan(u_0_full)));

X_use = X_contract_full;
X_use(rn_bin_use==0) = NaN;
pi_0_use = pi_0_full;
pi_0_use(rn_bin_use==0) = NaN;
u_0_use = u_0_full;
u_0_use(rn_bin_use==0)=NaN;

X_sims = [X_use(:,:,1);X_use(:,:,2)];
pi_0_sims = [pi_0_use(:,:,1);pi_0_use(:,:,2)];
u_0_sims = [u_0_use(:,:,1);u_0_use(:,:,2)];

%% Estimate true probabilities by conditioning on initial volatility
if estimate_true == 1
    
    % Setup for spot draws
    rng(seednum,"twister"); % rng seed
    S_sims = zeros(sim_T + 1,n); % each col. has one sim (burn-in plus true plus init)
    S_sims(1,:) = S_0;
    h_sims = zeros(sim_T + 1,n);
    h_sims(1,:) = h_0;
    z_sims = [zeros(1,n);randn(sim_T,n)];
    
    % RN volatility and price paths
    for t=2:sim_T+1            
        h_sims(t,:) = omega + beta*h_sims(t-1,:) + alpha*(z_sims(t-1,:)-gamma*sqrt(h_sims(t-1,:))).^2;
        S_sims(t,:) = exp(log(S_sims(t-1,:)) + rf + (mu-1/2)*h_sims(t,:) + sqrt(h_sims(t,:)).*z_sims(t,:));
    end
    
    % Realized return bins
    realized_ret = log(S_sims(exp_inds(2:end),:)./S_sims(exp_inds(1:end-1),:))-rf*T; % realized excess returns over option contracts
    max_return = 1; % min/max excess return capped at 100% for numerical integration
    bin_edges = [-max_return -0.2:0.05:0.2 max_return]; % discretize into 5-ppt ranges
    realized_bins = discretize(realized_ret,bin_edges);
    
    % Initial vol
    dt_inds_init = exp_inds(1:end-1);
    vol_0 = h_sims(dt_inds_init,:);

    % Other matrices as needed (following svcj file naming convention)
    realized_T = realized_bins;
    rn_bin_histogram = rn_bin_num_full;

    % Vectorize
    vol_full = reshape(vol_0,n_full*(exp_count-1),1);
    realized_full = reshape(realized_T,n_full*(exp_count-1),1);
    pi_0_star_binary_full = reshape(pi_0_full,n_full*(exp_count-1),2);
    rn_bin_full = reshape(rn_bin_histogram,n_full*(exp_count-1),3);
    [vol_cdf,vol_sorted] = ecdf(vol_full);
    pp = (0:0.01:1)'; % estimate probabilities across volatility state space, for 100 equally-spaced (by probability) ranges of volatility)
    vol_ends = NaN(size(pp));
    num_bins = max(max(rn_bin_full));
    realized_hist_by_vol = NaN(size(pp,1)-1,num_bins);
    pi_0_star_avg_by_vol = NaN(size(pp,1)-1,num_bins-1); % will be binary given format of saved rn probabilities above
    pi_0_star_sd_by_vol = NaN(size(pp,1)-1,num_bins-1);
    phi_by_vol = NaN(size(pp,1)-1,num_bins-1);
    binary_star_by_vol = pi_0_star_avg_by_vol;
    binary_true_by_vol = binary_star_by_vol;
    for ci = 1:numel(pp)
        vol_ends(ci) = vol_sorted(find(vol_cdf-pp(ci)<=0,1,'last'));
        pp(ci) = vol_cdf(find(vol_cdf-pp(ci)<=0,1,'last'));
        if ci>1
            vol_select = vol_full>=vol_ends(ci-1) & vol_full<vol_ends(ci);
            realized_use = realized_full(vol_select);
            realized_hist = histcounts(realized_use,(0.5:1:num_bins+0.5),'Normalization','probability');
            realized_hist_by_vol(ci-1,:) = realized_hist;
            pi_0_star_vol_use = reshape(pi_0_star_binary_full(vol_select,:),[],1);
            rn_bin_vol_use = reshape(rn_bin_full(vol_select,1:2),[],1);
            for bb = 1:num_bins-1
                binary_star_by_vol(ci-1,bb) = mean(pi_0_star_vol_use(rn_bin_vol_use==bb));
                pi_0_star_sd_by_vol(ci-1,bb) = std(pi_0_star_vol_use(rn_bin_vol_use==bb));
                binary_true_by_vol(ci-1,bb) = realized_hist_by_vol(ci-1,bb)/(realized_hist_by_vol(ci-1,bb)+realized_hist_by_vol(ci-1,bb+1));
                star_ratio = binary_star_by_vol(ci-1,bb)/(1-binary_star_by_vol(ci-1,bb));
                true_ratio = binary_true_by_vol(ci-1,bb)/(1-binary_true_by_vol(ci-1,bb));
                phi_by_vol(ci-1,bb) = star_ratio/true_ratio;
            end
        end
    end
end

%% Compute true phi by simulation and by bin
vol_pct = discretize(vol_0,vol_ends);
vol_pct_2 = [vol_pct;vol_pct];
calculated_phi_sims = NaN(size(X_sims));
rn_bin_sims = [rn_bin_pair(:,:,1);rn_bin_pair(:,:,2)];
X_all = [X_contract_full(:,:,1);X_contract_full(:,:,2)];
pi_0_all = [pi_0_full(:,:,1);pi_0_full(:,:,2)];
u_0_all = [u_0_full(:,:,1);u_0_full(:,:,2)];
max_bin = max(max(max(rn_bin_pair)));
X_by_bin = NaN(size(X_sims,1),size(X_sims,2),max_bin);
true_phi_bins = NaN(size(X_sims,1),size(X_sims,2),max_bin);
pi_0_by_bin = X_by_bin;
u_0_by_bin = X_by_bin;
for b=1:max_bin
    b_inds = rn_bin_sims==b;
    for j=1:n_full
        for i=1:(exp_count-1)*2
            if b_inds(i,j)==1
                calculated_phi_sims(i,j) = phi_by_vol(vol_pct_2(i,j),b);
                X_by_bin(i,j,b) = X_all(i,j);
                pi_0_by_bin(i,j,b) = pi_0_all(i,j);
                u_0_by_bin(i,j,b) = u_0_all(i,j);
                true_phi_bins(i,j,b) = calculated_phi_sims(i,j);
            end
            if isnan(X_sims(i,j))
                calculated_phi_sims(i,j) = NaN;
            end
        end
    end
end

%% Estimate phi and gamma using bounds
phi_sims = NaN(n,1);
phi_sims_alt = NaN(n,1);
true_phi_sims = NaN(n,1);
true_phi_bar_sims = NaN(n,1);
phi_by_bin = NaN(n,max_bin);
true_phi_by_bin = NaN(n,max_bin);
true_phi_bar_by_bin = NaN(n,max_bin);
fzopt = optimset('Display','off');
for j=1:n
    
    phi0 = 5;
    X_sim = X_sims(:,j);
    pi_0_sim = pi_0_sims(:,j);
    calculated_phi_sim = calculated_phi_sims(:,j);
    bound = @(phi) mean(X_sim-(1-(1./(phi+(1-phi)*pi_0_sim))).*pi_0_sim.^2,"omitnan");
    phi_sim = fzero(bound,phi0,fzopt);
    bound2 = @(phi) mean(-X_sim-(1-(1./(phi+(1-phi)*pi_0_sim))).*pi_0_sim.*(1-pi_0_sim),"omitnan");
    phi_sim2 = fzero(bound2,phi0,fzopt);
    phi_sims_alt(j) = max([phi_sim phi_sim2 1]);
    phi_sims(j) = max([phi_sim 1]);
    true_phi_sims(j) = mean(calculated_phi_sim,"omitnan");
    true_phi_bar_sims(j) = max(calculated_phi_sim);

    for b=1:max_bin
        X_sim_bin = X_by_bin(:,j,b);
        pi_0_sim_bin = pi_0_by_bin(:,j,b);
        phi_sim_bin = true_phi_bins(:,j,b);
        if isnan(mean(X_sim_bin,"omitnan"))
            continue;
        end
        if mean(X_sim_bin,"omitnan")>0
            boundb = @(phi) mean(X_sim_bin-(1-(1./(phi+(1-phi)*pi_0_sim_bin))).*pi_0_sim_bin.^2,"omitnan");
        else
            boundb = @(phi) mean(-X_sim_bin-(1-(1./(phi+(1-phi)*pi_0_sim_bin))).*pi_0_sim_bin.*(1-pi_0_sim_bin),"omitnan");
        end
        [phi_simb,~,flag1] = fzero(boundb,phi0,fzopt);
        phi_by_bin(j,b) = max([phi_simb 1]);
        true_phi_by_bin(j,b) = mean(phi_sim_bin(isfinite(phi_sim_bin)),"omitnan");
        true_phi_bar_by_bin(j,b) = max(phi_sim_bin(isfinite(phi_sim_bin)));
    end
end

%% Aggregate results

phi_overall = mean(phi_sims);
phi_overall_alt = mean(phi_sims_alt);
phi_bins = mean(phi_by_bin,1);
true_phi = mean(true_phi_sims);
true_phi_bar = mean(true_phi_bar_sims);
true_phi_bins = mean(true_phi_by_bin,1);
true_phi_bar_bins = mean(true_phi_bar_by_bin,1);

gamma_overall = (phi_overall-1)*20;
gamma_overall_alt = (phi_overall_alt-1)*20;
gamma_bins = (phi_bins-1)*20;
true_gamma_overall = (true_phi-1)*20;
true_gamma_bar = (true_phi_bar-1)*20;
true_gamma_bins = (true_phi_bins-1)*20;
true_gamma_bar_bins = (true_phi_bar_bins-1)*20;
gamma_true = mu;
true_rra1 = gamma_true;

%% Output table
X_bar = mean(mean(X_sims,1,"omitnan"));
u_0_bar = mean(mean(u_0_sims,1,"omitnan"));
X_u0_alt = mean(mean(X_sims,1,"omitnan")./mean(u_0_sims,1,"omitnan"));
X_over_u = round(X_bar/u_0_bar,2);
X_over_T = round(X_bar/T,4);
phi_bar_bound = round(phi_overall,1);
gamma_overall_bound = round(gamma_overall,1);
phi_bar_true = round(true_phi_bar,1);
gamma_bar_true = round(true_gamma_bar,1);
garch_table = table(X_over_u, X_over_T, phi_bar_bound,gamma_overall_bound,phi_bar_true,gamma_bar_true,'RowName',{'CHJ Model'});
disp(garch_table);
writetable(garch_table,'Output/table4_part1.txt','Delimiter','\t','WriteRowNames',true);
writetable(garch_table,'../FinalTables/table4_part1.txt','Delimiter','\t','WriteRowNames',true);

clearvars S_sims h_sims z_sims;
save("Output/Simulations/garch_output.mat","-v7.3");
