%% SVCJ model simulations (Table IV, part 2)
%
% Simulates the stochastic volatility with correlated jumps (SVCJ) 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_part2.txt, ../FinalTables/table4_part2.txt
%         Output/Simulations/svcj_sim_streams.mat, svcj_output.mat

clear
close all;
toolboxsubpath = 'HelperScripts/Toolbox';
addpath(genpath(toolboxsubpath));
addpath(genpath('Output'));
resimulate = 1; % set to 1 to rerun sims
recombine = 1; % set to 1 to reload all underlying sims data to combine and save
run_cpu = 1; % set to 0 if you're running with gpu; set to 1 if you need to run with cpu (much slower, so set N_full to some lower number than 10000 or else your sims will be running for months)
run_mex = 1*run_cpu; % if you have matlab coder installed and are running on cpu, runs faster if you create mex files. change to 0 if desired.
numThreads = 7; % I use 7 of my 8 logical cores for cpu parallelization; you can change if desired. for gpu, it's set to automatically use # of gpus below.

% Simulation parameters
N_full = 10000*(1-run_cpu) + 100*run_cpu; % simulations (I use 10000 for gpu, 100 for illustration with cpu)
n_tot = N_full/5; % sub-simulations for memory (set to 2000, equivalent to N_full/5 which also works if running fewer sims on cpu)
N = N_full/n_tot;
freq = 1; % MC steps of 1 day
dpy = 252; % trading days per year for simulations (note that Nicola's toolbox uses 365-day convention, but the underlying function calls will handle that for annualized inputs below)
dpm = 21; % per month
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 = 209; % picked randomly
rng(seednum,"twister"); % rng seed
T = 3*dpm; % 3-month contracts
T1 = T/dpy; % in years
exp_dates = zeros(1,sim_T + 1);
exp_dates(burn+1:T:end) = 1; % vector selecting expiration dates
exp_count = sum(exp_dates>0);
exp_inds = find(exp_dates>0);
rf = 0.05; % constant, following other simulations
d = 0.0; % no dividend yield

% Option pricing toolbox inputs
Nf = 500; % number of Fourier cosine approximation points (following toolbox recommendation)
Lf = 20; % truncation point for density computation in # of standard deviations (following toolbox recommendation)
distr = '2F_Gaussian'; % toolbox model (2F Gaussian encompasses a 1F double-jump model)
CP_use = 1; % will calculate RN CDF from digital call options (CP=1)
 
% Parameters under P [from Broadie Chernov Johannes 2007 Table I, SVCJ (row EJP, with exact parameters from that cited paper), annualized]
x_P = zeros(16,1); % preallocate parameter vector
x_P(1)   = -0.48;     % rho1 (parameter notation follows toolbox convention)
x_P(2)   = 0.0144;    % v1 (annualized as (sqrt(.5376*252)/100)^2)
x_P(3)   = 4.0320;    % k1 (252*0.0260)
x_P(4)   = 0.2;       % s1 (.079*252/100)
x_P(5)   =  0;    % rho2 (this model doesn't have a second vol factor)
x_P(6)   =  0;    % v2
x_P(7)   =  0;     % k2
x_P(8)   =  0;    % s2
x_P(9)   =  1.0080;   % c0 (jump intensity, .0066*252)
x_P(10)  =  0;     % c1
x_P(11)  =  0;      % c2
x_P(12)  = -0.0284;     % mu_y
x_P(13)  =  .0490;    % s_y
x_P(14)  = 0.0315;    % mu_v ((sqrt(1.48*252)/100)^2)
x_P(15)  = 0;      % rho_j (as in BCJ)
eta_s = 2.5; 
mu_y_bar_P = exp(x_P(12)+x_P(13)^2/2)-1;
mu_y_Q = -0.0501;
mu_y_bar_Q = exp(mu_y_Q+.0751^2/2)-1;
x_P(16) = x_P(9)*(mu_y_bar_P-mu_y_bar_Q); % equity premium (gamma_t in BCJ)

% Parameters under Q [also from BCJ: Table IV, last row (discussed in text as preferred model, with most reasonable parameter estimates)]
x_Q = x_P(1:15); % start from P vector (excluding equity premium) and modify as needed
x_Q(12)  = -0.0501;     % mu_y
x_Q(13)  =  .0751;    % s_y
x_Q(14)  = 0.0930;    % mu_v ((sqrt(3.71*252)/100)^2)

% Initialization parameters for start of burn-in period
S0 = 1; % spot price
sv1 = x_P(2); % variance at long-run value to start
sv2 = 0;
sv3 = 0; 

% Generate mex code if run_mex==1
if run_mex==1
    try % try an example call of the mex version of the digital option pricer -- will throw an error if either the mex file doesn't exist or is incompatible with the local machine
        S_use = 1.2; % setting some arbitrary parameters
        K_use = 1.25;
        ytm_j = 0.25;
        sv1_use = sv1;
        [~,rn_out,~,~,~] = digital_option_price_cg_mex(distr,S_use,K_use,x_Q,ytm_j,CP_use,sv1_use,sv2,sv3,rf,d,Nf,Lf);
        use_mex = 1;
    catch
        cd(toolboxsubpath);
        try
            build_digital_mex(distr,x_Q);
            use_mex=1;
        catch
            use_mex=0;
        end
        cd ../..;
        addpath(genpath(toolboxsubpath));
    end
end

%% Simulation runs
if resimulate==1
    
    if run_cpu==1 && isempty(gcp('nocreate')) % try to parallelize
        % Thread-based pool avoids copying big arrays to worker processes
        parpool("Processes",numThreads);
    elseif run_cpu==0 && isempty(gcp('nocreate')) % for gpu case, only parallelize with as many threads as gpus 
        parpool("Threads",gpuDeviceCount);
    end

    % Sub-simulations
    for n_i = 1:n_tot
        tic;
        disp(n_i);

        % Setup for simulation draws, following toolbox notational conventions
        S  = zeros(N,sim_T+1); % now each row has one sim (burn-in plus true plus init)
        V1 = zeros(N,sim_T+1);
        V2 = zeros(N,sim_T+1);
        S(:,1) = S0;
        V1(:,1) = sv1;
        V2(:,1) = sv2;
        dt = 1/(dpy*freq); % annualized time step
        x = x_P; % relabel to plug into Nicola's code (simulate under P)
        rho1 = x(1); v1 = x(2); k1 = x(3); s1   = x(4); rho2 = x(5); v2 = x(6); k2 = x(7); s2 = x(8); c0 = x(9); c1 = x(10); c2 = x(11); muy = x(12); sy = x(13); mv = x(14); rhoj = x(15); gam = x(16);
        comp =  (exp(muy+0.5*sy^2)/(1-rhoj*mv)) - 1.0;   % jump compensator E[exp(-1*X)] 
    
        % Volatility and price path draws, using toolbox code
        for tt=2:sim_T+1
            
            % Gaussian shocks for this day
            eps1 = randn(N,1);
            eps2 = randn(N,1);
            eps3 = randn(N,1);
            eps4 = randn(N,1);
            
            % Jump draws
            c = c0; % don't need c1 or c2 given constant jump intensity
            lmd = c*dt*ones(N,1); % constant jump intensity
            Pj = zeros(N,1); % jump vectors
            Jv = zeros(N,1);
            Jx = zeros(N,1);
            for j=1:N
                Pj(j,1) = getBinomial(1,min(lmd(j),1));
                if Pj(j,1) == 1
                    Jv(j,1) = exprnd(mv,1,1);
                    Jx(j,1) = normrnd(muy+rhoj*Jv(j,1),sy,1,1);                       
                end
            end
            
            % Realized log return, volatility, and price steps
            eqp = eta_s.*V1(:,tt-1)+gam;
            r = ((rf-d+eqp)-0.5*(V1(:,tt-1)) - c*comp)*dt + sqrt(V1(:,tt-1)*dt).*eps1(:,1) + Jx;
            V1(:,tt) = abs(V1(:,tt-1) + k1.* (v1-V1(:,tt-1))*dt + s1.*sqrt(V1(:,tt-1).*dt).*(rho1*eps1(:,1)+sqrt(1-rho1^2)*eps2(:,1))) + Jv;
            S(:,tt) = S(:,tt-1).*exp(r);
        end

        clearvars V2 Jv Jx Pj eps1 eps2 eps3 eps4 lmd;

        % Setup for option price calcs
        realized_ret = log(S(:,exp_inds(2:end))./S(:,exp_inds(1:end-1))) - rf*T1; % realized excess returns over option contracts
        m_max = 10; % min/max excess return (10 standard deviations) for numerical integration
        max_return = m_max*sqrt(sv1)*sqrt(T1); % min/max excess return cap
        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(:,exp_inds(1:end-1));
        start_array = repmat(start_level,[1 1 4]);
        strikes_pre = start_array.*exp(rf*T1+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(N,exp_count-1,4,T); % dims: 1 for simulation number, 2 for contract #, 3 for strike, 4 for trading day within contract
        dtm = repmat(shiftdim(T:-1:1,-2),[N,exp_count-1,4,1]);
        prices = NaN(size(dtm));
        vols = NaN(size(dtm));
        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;
        exp_inds_use = exp_inds(2:end);

        % Run through by dtm
        parfor j=1:T
            dtm_j = T-j+1;
            ytm_j = dtm_j/dpy; % years to maturity
            dt_inds = exp_inds_use-dtm_j;
            prices_j = S(:,dt_inds);
            vols_j = V1(:,dt_inds);
            price_slice = repmat(prices_j,[1 1 4]);
            strike_slice = strikes(:,:,:,j);
            vol_slice = repmat(vols_j,[1 1 4]);
            
            % If running on GPU
            if run_cpu==0
                % Set up for GPU call
                prices_gpu = gpuArray(price_slice);
                strikes_gpu = gpuArray(strike_slice);
                vols_gpu = gpuArray(vol_slice);
        
                % Need scalar or dimension-matching inputs for GPU arrayfun
                x = x_Q; % copied from relabeling above for P parameters
                rho1 = x(1); v1 = x(2); k1 = x(3); s1 = x(4); rho2 = x(5); v2 = x(6); k2 = x(7); s2 = x(8); c0 = x(9); c1 = x(10); c2 = x(11); muy = x(12); sy = x(13); mv = x(14); rhoj = x(15);
                
                % Compute RN CDF for each simulation and expiration date
                rn_out = arrayfun(@digital_option_price_gpu,prices_gpu,strikes_gpu,ytm_j,CP_use,vols_gpu,sv2,sv3,rf,d,Nf,Lf,rho1,v1,k1,s1,rho2,v2,k2,s2,c0,c1,c2,muy,sy,mv,rhoj);
                rn_mat = gather(rn_out);
                rn_cdf(:,:,:,j) = rn_mat;

            else % if running on CPU
                
                if run_mex==1 && use_mex==1 && exist('digital_option_price_cg_mex','file')
                    pricer = @digital_option_price_cg_mex;
                else
                    pricer = @digital_option_price;
                end

                rn_cdf_slice = NaN(N,exp_count-1,4,1);
                for n=1:N
                    for exp_j=1:exp_count-1
                        sv1_use = vols_j(n,exp_j);
                        S_use = squeeze(price_slice(n,exp_j,:));
                        K_use = squeeze(strike_slice(n,exp_j,:));
                        [~,rn_out,~,~,~] = pricer(distr,S_use,K_use,x_Q,ytm_j,CP_use,sv1_use,sv2,sv3,rf,d,Nf,Lf);
                        rn_cdf_slice(n,exp_j,:,1) = rn_out;
                    end
                end
                rn_cdf(:,:,:,j) = rn_cdf_slice;
            end
        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),[N exp_count-1 1 1]); % for last day
        rn_probs = cat(4,rn_probs,realization);
        dt_inds_init = exp_inds_use-T;
        vols_j = V1(:,dt_inds_init); % use for initial vol below
        if run_cpu==0
            file_init = "Output/Simulations/svcj_sims_use_gpu_%d.mat";
        else
            file_init = "Output/Simulations/svcj_sims_use_%d.mat";
        end
        file_nm = sprintf(file_init,n_i);
	save(file_nm,'-v7.3');
	toc;
    end
end

%% Load and combine simulations
if recombine==1
    if run_cpu==0
        file_init = "Output/Simulations/svcj_sims_use_gpu_%d.mat";
    else
        file_init = "Output/Simulations/svcj_sims_use_%d.mat";
    end
    pi_0_star_histogram = NaN(N_full,exp_count-1,3);
    rn_bin_histogram = NaN(N_full,exp_count-1,3);
    vol_0 = NaN(N_full,exp_count-1);
    realized_T = NaN(N_full,exp_count-1);
    X_contract_full = NaN(N_full,exp_count-1,2);
    X_use_full = NaN(N_full,exp_count-1,2);
    pi_0_use_full = NaN(N_full,exp_count-1,2);
    pi_0_full = NaN(N_full,exp_count-1,2);
    u_0_full = NaN(N_full,exp_count-1,2);
    u_0_use_full = NaN(N_full,exp_count-1,2);
    rn_bin_pair_full = NaN(N_full,exp_count-1,2);
    rn_bin_use_full = NaN(N_full,exp_count-1,2);
    for n_i=1:n_tot
        
        % Load and clean
        file_nm = sprintf(file_init,n_i);
        load(file_nm,"rn_probs","rn_bin_num","vols_j","realized_bins");
        rn_probs(rn_probs<0)=0;
        rn_probs(rn_probs>1)=1;
        n_inds = (n_i-1)*N+1:n_i*N;
    
        % Store RN probabilities, realizations, starting volatility to estimate true phi
        pi_0_star_histogram(n_inds,:,:) = rn_probs(:,:,:,1);
        rn_bin_histogram(n_inds,:,:) = rn_bin_num;
        vol_0(n_inds,:) = vols_j;
        realized_T(n_inds,:) = realized_bins;
    
        % Binarize beliefs and create main binary statistics
        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);
        rn_bin_pair = rn_bin_num(:,:,1:2);
        rn_bin_pair(isnan(u_0)) = NaN;
        rn_bin_use = ~isnan(rn_bin_pair) & rn_bin_pair>=2 & rn_bin_pair<=8;
        X_use = X_contract;
        X_use(rn_bin_use==0) = NaN;
        pi_0_use = pi_0;
        pi_0_use(rn_bin_use==0) = NaN;
        u_0_use = u_0;
        u_0_use(rn_bin_use==0) = NaN;
    
        % Combine binarized statistics
        X_contract_full(n_inds,:,:) = X_contract;
        X_use_full(n_inds,:,:) = X_use;
        pi_0_full(n_inds,:,:) = pi_0;
        pi_0_use_full(n_inds,:,:) = pi_0_use;
        u_0_full(n_inds,:,:) = u_0;
        u_0_use_full(n_inds,:,:) = u_0_use;
        rn_bin_pair_full(n_inds,:,:) = rn_bin_pair;
        rn_bin_use_full(n_inds,:,:) = rn_bin_use;
    
    end
    
    % Clean and save
    X_contract = X_contract_full;
    X_use = X_use_full;
    pi_0 = pi_0_full;
    pi_0_use = pi_0_use_full;
    u_0 = u_0_full;
    u_0_use = u_0_use_full;
    rn_bin_pair = rn_bin_pair_full;
    rn_bin_use = rn_bin_use_full;
    clearvars X_contract_full X_use_full pi_0_full pi_0_use_full u_0_full u_0_use_full rn_bin_pair_full rn_bin_use_full
    save("Output/Simulations/svcj_sim_streams.mat","X_contract","X_use","pi_0","pi_0_use","u_0","u_0_use","rn_bin_pair","rn_bin_use","pi_0_star_histogram","rn_bin_histogram","vol_0","realized_T","-v7.3");
else
    load("Output/Simulations/svcj_sim_streams.mat");
end

%% Estimate true probabilities using all sims
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_full = reshape(pi_0_star_histogram,N_full*(exp_count-1),3);
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);
pi_0_star_sd_by_vol = NaN(size(pp,1)-1,num_bins);
phi_by_vol = NaN(size(pp,1)-1,num_bins-1);
binary_star_by_vol = NaN(size(pp,1)-1,num_bins-1);
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_full(vol_select,:),[],1);
        rn_bin_vol_use = reshape(rn_bin_full(vol_select,:),[],1);
        for bb = 1:num_bins
            pi_0_star_avg_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));
            if bb>1
                binary_star_by_vol(ci-1,bb-1) = pi_0_star_avg_by_vol(ci-1,bb-1)/(pi_0_star_avg_by_vol(ci-1,bb-1)+pi_0_star_avg_by_vol(ci-1,bb));
                binary_true_by_vol(ci-1,bb-1) = realized_hist_by_vol(ci-1,bb-1)/(realized_hist_by_vol(ci-1,bb-1)+realized_hist_by_vol(ci-1,bb));
                star_ratio = binary_star_by_vol(ci-1,bb-1)/(1-binary_star_by_vol(ci-1,bb-1));
                true_ratio = binary_true_by_vol(ci-1,bb-1)/(1-binary_true_by_vol(ci-1,bb-1));
                phi_by_vol(ci-1,bb-1) = star_ratio/true_ratio;
            end
        end
    end
end

%% Compute true phi by simulation and by bin
X_sims = [X_use(:,:,1) X_use(:,:,2)];
vol_pct = discretize(vol_0,vol_ends);
vol_pct_2 = [vol_pct vol_pct];
pi_0_sims = [pi_0_use(:,:,1) pi_0_use(:,:,2)];
calculated_phi_sims = NaN(size(X_sims));
u_0_sims = [u_0_use(:,:,1) u_0_use(:,:,2)];
rn_bin_sims = [rn_bin_pair(:,:,1) rn_bin_pair(:,:,2)];
X_all = [X_contract(:,:,1) X_contract(:,:,2)];
pi_0_all = [pi_0(:,:,1) pi_0(:,:,2)];
u_0_all = [u_0(:,:,1) u_0(:,:,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 i=1:N_full
        for j=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_full,1);
phi_sims_alt = NaN(N_full,1);
true_phi_sims = NaN(N_full,1);
true_phi_bar_sims = NaN(N_full,1);
phi_by_bin = NaN(N_full,max_bin);
true_phi_by_bin = NaN(N_full,max_bin);
for j=1:N_full
    
    phi_0 = 5;
    X_sim = X_sims(j,:);
    pi_0_sim = pi_0_sims(j,:);
    calculated_phi_sim = calculated_phi_sims(j,:);
    bound1 = @(phi) mean(X_sim-(1-(1./(phi+(1-phi)*pi_0_sim))).*pi_0_sim.^2,"omitnan");
    phi_sim1 = fzero(bound1,phi_0);
    bound2 = @(phi) mean(-X_sim-(1-(1./(phi+(1-phi)*pi_0_sim))).*pi_0_sim.*(1-pi_0_sim),"omitnan");
    opts2 = optimset('display', 'none');
    phi_sim2 = fzero(bound2,phi_0,opts2);
    phi_sims_alt(j) = max([phi_sim1 phi_sim2 1]);
    phi_sims(j) = max([phi_sim1 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,phi_0);
        phi_by_bin(j,b) = max([phi_simb 1]);
        true_phi_by_bin(j,b) = mean(phi_sim_bin,"omitnan");
    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);

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_rra1 = (x_P(12)-mu_y_Q)/x_P(13)^2; % from BCJ, citing Bates (1988) and assuming covariance of 1 with stock-market jumps and wealth jumps

%% Output table
X_bar = mean(mean(X_sims,2,"omitnan"));
u_0_bar = mean(mean(u_0_sims,2,"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);
svcj_table = table(X_over_u, X_over_T, phi_bar_bound,gamma_overall_bound,phi_bar_true,gamma_bar_true,'RowName',{'SVCJ Model'});
disp(svcj_table);
writetable(svcj_table,'Output/table4_part2.txt','Delimiter','\t','WriteRowNames',true);
writetable(svcj_table,'../FinalTables/table4_part2.txt','Delimiter','\t','WriteRowNames',true);
save("Output/Simulations/svcj_output.mat","-v7.3");