%%%%%%%%%%%%%% NOTE: TO RUN THIS PROGRAM (if error occurs) %%%%%%%%%%%%%%%%%%%
% - When running Matlab for creating basis functions, 
% need to evoke the following commands after installing COMPECON:

tic
cepath='C:\Users\domnisc1\Dropbox\structural econometrics course not shared\part 2\COMPECON\';
path([cepath 'cetools;' cepath 'cedemos'],path);


clear all
load 'C:\Users\domnisc1\Dropbox\structural econometrics course not shared\part 2\Banking_RFS15\data_files\ccpdata_pvar_KLW.mat'

aftersimulationfix = 1;
%{ 
  House_index              48245x1             385960  double              
  House_index_next         48245x1             385960  double              
  Senate_index             48245x1             385960  double              
  Senate_index_next        48245x1             385960  double              
  V                           11x11               968  double              
  assGrowth_1y_lag1        48245x1             385960  double              
  assGrowth_1y_lag2        48245x1             385960  double              
  assGrowth_1y_lag3        48245x1             385960  double              
  assGrowth_1y_lag4        48245x1             385960  double              
  asset_growth_1y          48245x1             385960  double              
  asset_growth_2y          48245x1             385960  double              
  assets                   48245x1             385960  double              
  assets_lag1              48245x1             385960  double              
  assets_lag2              48245x1             385960  double              
  assets_lag3              48245x1             385960  double              
  assets_lag4              48245x1             385960  double              
  equity                   48245x1             385960  double              
  equity_a                 48245x1             385960  double              
  equity_a_lag1            48245x1             385960  double              
  equity_a_lag2            48245x1             385960  double              
  equity_a_lag3            48245x1             385960  double              
  equity_a_lag4            48245x1             385960  double              
  estimated_cost           48245x1             385960  double              
  failed                   48245x1             385960  double              
  first_core_dep_a         48245x1             385960  double              
  log_assets               48245x1             385960  double              
  log_assets_lag1          48245x1             385960  double              
  log_assets_lag2          48245x1             385960  double              
  log_assets_lag3          48245x1             385960  double              
  log_assets_lag4          48245x1             385960  double              
  log_county_deposits      48245x1             385960  double              
  net_income               48245x1             385960  double              
  net_income_annual        48245x1             385960  double              
  nonperf_a_lag1           48245x1             385960  double              
  nonperf_a_lag2           48245x1             385960  double              
  nonperf_a_lag3           48245x1             385960  double              
  nonperf_a_lag4           48245x1             385960  double              
  nonperf_loans            48245x1             385960  double              
  nonperf_loans_a          48245x1             385960  double              
  realEst_own_a_lag1       48245x1             385960  double              
  realEst_own_a_lag2       48245x1             385960  double              
  realEst_own_a_lag3       48245x1             385960  double              
  realEst_own_a_lag4       48245x1             385960  double              
  real_estate_owned        48245x1             385960  double              
  real_estate_owned_a      48245x1             385960  double              
  roa_annual               48245x1             385960  double              
  roa_annual_lag1          48245x1             385960  double              
  roa_annual_lag2          48245x1             385960  double              
  roa_annual_lag3          48245x1             385960  double              
  roa_annual_lag4          48245x1             385960  double              
  roa_qtr                  48245x1             385960  double              
  rssd_id                  48245x1             385960  double              
  state                    48245x1             385960  double              
  unemp_nat                48245x1             385960  double              
  unemp_state              48245x1             385960  double              
  unemp_state_lag1         48245x1             385960  double              
  unemp_state_lag2         48245x1             385960  double              
  unemp_state_lag3         48245x1             385960  double              
  unemp_state_lag4         48245x1             385960  double              
  us_equity_a              48245x1             385960  double              
  us_equity_a_lag1         48245x1             385960  double              
  yearq                    48245x1             385960  double    
%}


ALLX = [rssd_id, yearq, failed, state, ... %1 to 4
  equity_a, equity_a_lag1, equity_a_lag2, equity_a_lag3, ... % 5 to 8
  nonperf_loans_a, nonperf_a_lag1, nonperf_a_lag2, nonperf_a_lag3, ... % 9 to 12
  log_assets, log_assets_lag1, log_assets_lag2, log_assets_lag3, ... % 13 to 16
  real_estate_owned_a, realEst_own_a_lag1, realEst_own_a_lag2, realEst_own_a_lag3, ...% 17 to 20
  roa_annual, roa_annual_lag1, roa_annual_lag2, roa_annual_lag3, ...% 21 to 24
  unemp_state, unemp_state_lag1, unemp_state_lag2, unemp_state_lag3,... % 25 to 28
  asset_growth_1y, assGrowth_1y_lag1, assGrowth_1y_lag2, assGrowth_1y_lag3, ... % 29 to 32
  House_index, Senate_index, House_index_next, Senate_index_next... % 33, 34, 35, 36
  assets, assets_lag1, assets_lag2, assets_lag3, ...; % 37 to 40
  first_core_dep_a,log_county_deposits,us_equity_a,us_equity_a_lag1]; % 41, 42, 43, 44

J = sum(isnan(ALLX(:,1:end-4)) + isinf(ALLX(:,1:end-4)),2); % note that you don't want to drop the observation if the last four variables are not available (this is ver.26 not ver.27)
J = logical(J > 0);
JJ = logical(net_income_annual >= 0); % no case
sum(JJ)
JJJ = logical(equity_a > 0.3); % no case
sum(JJJ)
JJJJ = logical(yearq > (1992 - 1960)*4 + (4-1) | yearq < (1985 - 1960)*4 + (4-1)); %only use between (1985Q4~1992Q4) 
% for transitions: "keep if inrange(yearq, yq(1985,4),yq(1992,4))"
% for cost function estimation: from 1985.3 to 1992.4
% after that, data range is 1985.4 to 1992.4
J = logical(J > 0 | JJ > 0  | JJJ > 0 | JJJJ > 0);
sum(J)

Jtemp = logical(J == 1 & ALLX(:,3) == 1); % loosing 190 closed banks since they don't have lags of varibles.
ALLXtemp = ALLX(Jtemp,:);

ALLX(J,:) = [];
ALLX = sortrows(ALLX,[1 2]);
BankID = ALLX(:,1);
yearq = ALLX(:,2);
Bclose = ALLX(:,3);
equity_a_all = ALLX(:,5:8);
npf_a_all = ALLX(:,9:12);
logasset_all = ALLX(:,13:16);
realEst_own_a_all = ALLX(:,17:20);
roa_annual_all = ALLX(:,21:24);
unemp_all = ALLX(:,25:28);%why do we need lags? they used in transition of other bank variables.
assGrowth_1y_all = ALLX(:,29:32);
House = ALLX(:,33);
Senate = ALLX(:,34); 
House_next = ALLX(:,35);
Senate_next = ALLX(:,36); 
asset_all = ALLX(:,37:40);
core_dep_a = ALLX(:,41);
logcounty_dep = ALLX(:,42);
us_equity_a_all = ALLX(:,43:44);

Jclose = logical(Bclose == 1);
Jopen = logical(Bclose == 0);

num_obs = size(ALLX,1);

beta_bankvar = xlsread('C:\Users\domnisc1\Dropbox\structural econometrics course not shared\part 2\Banking_RFS15\data_files\coeff_resid_KLW.xlsx','trans_beta');%note that matlab does not import Strings. So the first col. in the excel file is not existing.
beta_logasset = beta_bankvar(:,1);
J = isnan(beta_logasset);
beta_logasset(J,:) = [];
beta_equity_a = beta_bankvar(:,2);
J = isnan(beta_equity_a);
beta_equity_a(J,:) = [];
beta_npf_a = beta_bankvar(:,3);
J = isnan(beta_npf_a);
beta_npf_a(J,:) = [];
beta_realest_a = beta_bankvar(:,4);
J = isnan(beta_realest_a);
beta_realest_a(J,:) = [];
beta_roa = beta_bankvar(:,5);
J = isnan(beta_roa);
beta_roa(J,:) = [];
beta_grow_a = beta_bankvar(:,6);
J = isnan(beta_grow_a);
beta_grow_a(J,:) = [];
beta_asset = beta_bankvar(:,7);
J = isnan(beta_asset);
beta_asset(J,:) = [];

%{
L.log_assets
L2.log_assets
L3.log_assets
L4.log_assets
L.equity_a
L2.equity_a
L3.equity_a
L4.equity_a
L.nonperf_loans_a
L2.nonperf_loans_a
L3.nonperf_loans_a
L4.nonperf_loans_a
L.real_estate_owned_a
L2.real_estate_owned_a
L3.real_estate_owned_a
L4.real_estate_owned_a
L.roa_annual
L2.roa_annual
L3.roa_annual
L4.roa_annual
L.asset_growth_1y
L2.asset_growth_1y
L3.asset_growth_1y
L4.asset_growth_1y
L.unemp_state
L2.unemp_state
L3.unemp_state
L4.unemp_state
%}
Xfortran = [logasset_all,equity_a_all,npf_a_all,...%remember it is for tomorrow's, so we need to use current as the first lag.
realEst_own_a_all,roa_annual_all,assGrowth_1y_all,unemp_all,ones(num_obs,1)];

Elogasset_next = Xfortran*beta_logasset;
Eequity_a_next = Xfortran*beta_equity_a;
Enpf_a_next = Xfortran*beta_npf_a;
Erealest_a_next = Xfortran*beta_realest_a;
Eroa_next = Xfortran*beta_roa;
Egrow_a_next = Xfortran*beta_grow_a;
Easset_next = [asset_all,ones(num_obs,1)]*beta_asset;


bankvar_residuals = xlsread('C:\Users\domnisc1\Dropbox\structural econometrics course not shared\part 2\Banking_RFS15\data_files\coeff_resid_KLW.xlsx','trans_res');
residuals_logasset = bankvar_residuals(:,1+2); % _2 is to take into account the first two columns of bankid and yearq.
J = isnan(residuals_logasset);
residuals_logasset(J,:) = [];
residuals_equity_a = bankvar_residuals(:,2+2);
J = isnan(residuals_equity_a);
residuals_equity_a(J,:) = [];
residuals_npf_a = bankvar_residuals(:,3+2);
J = isnan(residuals_npf_a);
residuals_npf_a(J,:) = [];
residuals_realest_a = bankvar_residuals(:,4+2);
J = isnan(residuals_realest_a);
residuals_realest_a(J,:) = [];
residuals_roa = bankvar_residuals(:,5+2);
J = isnan(residuals_roa);
residuals_roa(J,:) = [];
residuals_grow_a = bankvar_residuals(:,6+2);
J = isnan(residuals_grow_a);
residuals_grow_a(J,:) = [];
residuals_asset = bankvar_residuals(:,7+2);
J = isnan(residuals_asset);
residuals_asset(J,:) = [];

beta_macrovar = xlsread('C:\Users\domnisc1\Dropbox\structural econometrics course not shared\part 2\Banking_RFS15\data_files\coeff_resid_KLW.xlsx','unemp_beta_AR4');
beta_unemp = beta_macrovar(:,1);
J = isnan(beta_unemp);
beta_unemp(J,:) = [];

Eunemp_next = beta_unemp(1)*unemp_all(:,1) + beta_unemp(2)*unemp_all(:,2) + beta_unemp(3)*unemp_all(:,3) + beta_unemp(4)*unemp_all(:,4) + beta_unemp(5)*ones(size(unemp_all(:,1)));

macrovar_residuals = xlsread('C:\Users\domnisc1\Dropbox\structural econometrics course not shared\part 2\Banking_RFS15\data_files\coeff_resid_KLW.xlsx','unemp_res_AR4');
residuals_unemp = macrovar_residuals(:,2); % the first column is yearq
J = isnan(residuals_unemp);
residuals_unemp(J,:) = [];

% ADD: To avoid extreme values of future probabilities or other future
% variables, the residuals are drawn from 0.25 percentile to 99.75
% percentile.

pctl_limit = prctile(residuals_logasset, [0.25 99.75]);
J1 = logical(residuals_logasset > pctl_limit(1) & residuals_logasset < pctl_limit(2));
pctl_limit = prctile(residuals_equity_a, [0.25 99.75]);
J2 = logical(residuals_equity_a > pctl_limit(1) & residuals_equity_a < pctl_limit(2));
pctl_limit = prctile(residuals_npf_a, [0.25 99.75]);
J3 = logical(residuals_npf_a > pctl_limit(1) & residuals_npf_a < pctl_limit(2));
pctl_limit = prctile(residuals_realest_a, [0.25 99.75]);
J4 = logical(residuals_realest_a > pctl_limit(1) & residuals_realest_a < pctl_limit(2));
pctl_limit = prctile(residuals_roa, [0.25 99.75]);
J5 = logical(residuals_roa > pctl_limit(1) & residuals_roa < pctl_limit(2));
pctl_limit = prctile(residuals_grow_a, [0.25 99.75]);
J6 = logical(residuals_grow_a > pctl_limit(1) & residuals_grow_a < pctl_limit(2));
pctl_limit = prctile(residuals_asset, [0.25 99.75]);
J7 = logical(residuals_asset > pctl_limit(1) & residuals_asset < pctl_limit(2));
pctl_limit = prctile(residuals_unemp, [0.25 99.75]);
J8 = logical(residuals_unemp > pctl_limit(1) & residuals_unemp < pctl_limit(2));
[f1,ee1] = ecdf(residuals_logasset(J1,:)); 
[f2,ee2] = ecdf(residuals_equity_a(J2,:));  
[f3,ee3] = ecdf(residuals_npf_a(J3,:));
[f4,ee4] = ecdf(residuals_realest_a(J4,:));
[f5,ee5] = ecdf(residuals_roa(J5,:));
[f6,ee6] = ecdf(residuals_grow_a(J6,:));
[f7,ee7] = ecdf(residuals_asset(J7,:));
[f8,ee8] = ecdf(residuals_unemp(J8,:));

beta_mc_sigma = xlsread('C:\Users\domnisc1\Dropbox\structural econometrics course not shared\part 2\Banking_RFS15\data_files\coeff_resid_KLW.xlsx','cost_beta_v26');
%{ 
Monetary Cost by v26.
model	
log_assets	-.0307782
equity_a	-.9750841
nonperf_loans_a	1.614266
real_estate_owned_a	1.377229
roa_annual	-.9223549
unemp_state	.0138776
_cons	.0758327
sigma	
_cons	.1605512
%}
beta_mc = beta_mc_sigma(1:end-1,:); sigma_mc = beta_mc_sigma(end);

Xmc_now = [logasset_all(:,1), equity_a_all(:,1), npf_a_all(:,1), realEst_own_a_all(:,1), roa_annual_all(:,1), unemp_all(:,1), ones(num_obs,1)];
mc_now_temp = Xmc_now*beta_mc;
mc_now_temp = normcdf(mc_now_temp/sigma_mc,0,1).*mc_now_temp+sigma_mc*normpdf(mc_now_temp/sigma_mc,0,1); % See Wooldridge p674 (Eq.17.14)
Emc_now = mc_now_temp.*asset_all(:,1);

AA = importdata('C:\Users\domnisc1\Dropbox\structural econometrics course not shared\part 2\Banking_RFS15\output_files\logit_beta_KLW.txt'); 
plogit_beta = AA.data; % note: it is from Stata=> the constant is the last one.
%{
    %%% If the model with trend is used:
AA = importdata('C:/Users/LoweryR/Documents/My Dropbox/UND_research_Banking/Final_files/output_files/logit_trend_beta_KLW.txt');
plogit_trend_beta = AA.data; % note: it is from Stata=> the constant is the last one.
%}
num_lags = 3; % for each observation, create basis functions for l1,l2 and l3 (current one will be made below while simulating, hence like in prob_param_v25, the true value is 4.)
num_basis = 4; % when 5 is used, it converges in 29 iterations but has some backed ups while using 4 knots it iterates in 15 iterations.
num_basis_a = num_basis-1; % we exclude the last columns to avoid multicolinearity (remember: sum(basis,2) == 1)
equity_a_all_basis_temp = NaN(length(equity_a_all(:,1)),num_basis_a*num_lags);
npf_a_all_basis_temp = NaN(length(equity_a_all(:,1)),num_basis_a*num_lags);
logasset_all_basis_temp = NaN(length(equity_a_all(:,1)),num_basis_a*num_lags);
realEst_own_a_all_basis_temp = NaN(length(equity_a_all(:,1)),num_basis_a*num_lags);
roa_annual_all_basis_temp = NaN(length(equity_a_all(:,1)),num_basis_a*num_lags);
assGrowth_1y_all_basis_temp = NaN(length(equity_a_all(:,1)),num_basis_a*num_lags);

% NOTE! we use min/max out of ALL samples. The knot locations are consistent across the lags.
min_equity_a = min(min(equity_a_all)); max_equity_a = max(max(equity_a_all));
min_npf_a = min(min(npf_a_all)); max_npf_a = max(max(npf_a_all));
min_logasset = min(min(logasset_all)); max_logasset = max(max(logasset_all));
min_realEst_own_a = min(min(realEst_own_a_all)); max_realEst_own_a = max(max(realEst_own_a_all));
min_roa_annual = min(min(roa_annual_all)); max_roa_annual = max(max(roa_annual_all));
min_assGrowth_1y = min(min(assGrowth_1y_all)); max_assGrowth_1y = max(max(assGrowth_1y_all));

for i = 1:num_lags % the last lag is not used, since we will have "next" as the first element.
fspace = fundefn('spli',num_basis, min_equity_a,max_equity_a);
equity_a_basis = funbas(fspace,equity_a_all(:,i));
equity_a_basis = full(equity_a_basis);
equity_a_all_basis_temp(:,num_basis_a*(i-1)+1:num_basis_a*(i-1)+num_basis_a) = equity_a_basis(:,1:end-1);
% basis functions for npf_a;
fspace = fundefn('spli',num_basis, min_npf_a,max_npf_a);
npf_a_basis = funbas(fspace,npf_a_all(:,i));
npf_a_basis = full(npf_a_basis);
npf_a_all_basis_temp(:,num_basis_a*(i-1)+1:num_basis_a*(i-1)+num_basis_a) = npf_a_basis(:,1:end-1);
% basis functions for logasset;
fspace = fundefn('spli',num_basis, min_logasset,max_logasset);
logasset_basis = funbas(fspace,logasset_all(:,i));
logasset_basis = full(logasset_basis);
logasset_all_basis_temp(:,num_basis_a*(i-1)+1:num_basis_a*(i-1)+num_basis_a) = logasset_basis(:,1:end-1);
% basis functions for realEst_own_a;
fspace = fundefn('spli',num_basis, min_realEst_own_a,max_realEst_own_a);
realEst_own_a_basis = funbas(fspace,realEst_own_a_all(:,i));
realEst_own_a_basis = full(realEst_own_a_basis);
realEst_own_a_all_basis_temp(:,num_basis_a*(i-1)+1:num_basis_a*(i-1)+num_basis_a) = realEst_own_a_basis(:,1:end-1);
% basis functions for roa_annual;
fspace = fundefn('spli',num_basis, min_roa_annual,max_roa_annual);
roa_annual_basis = funbas(fspace,roa_annual_all(:,i));
roa_annual_basis = full(roa_annual_basis);
roa_annual_all_basis_temp(:,num_basis_a*(i-1)+1:num_basis_a*(i-1)+num_basis_a) = roa_annual_basis(:,1:end-1);
% basis functions for assGrowth_1y;
fspace = fundefn('spli',num_basis, min_assGrowth_1y,max_assGrowth_1y);
assGrowth_1y_basis = funbas(fspace,assGrowth_1y_all(:,i));
assGrowth_1y_basis = full(assGrowth_1y_basis);
assGrowth_1y_all_basis_temp(:,num_basis_a*(i-1)+1:num_basis_a*(i-1)+num_basis_a) = assGrowth_1y_basis(:,1:end-1);
end

Eln_prob_Bclose = NaN(num_obs,1);
Eprob_Bclose = NaN(num_obs,1);
Eln_prob_trend_Bclose = NaN(num_obs,1);
Eprob_trend_Bclose = NaN(num_obs,1);
Emc_next = NaN(num_obs,1);
% Emc_next_pos_cost = NaN(num_obs,1);
Easset_next2 = NaN(num_obs,1);
Elogasset_next2 = NaN(num_obs,1);
Enpf_a_next2 = NaN(num_obs,1);
Eroa_next2 = NaN(num_obs,1);
Eunemp_state_next = NaN(num_obs,1);
% num_pos_cost = NaN(num_obs,1);
num_v = 8;
num_sim = 5000;

load 'C:\Users\domnisc1\Dropbox\structural econometrics course not shared\part 2\Banking_RFS15\data_files\rng_KLW.mat';




tic
for i = 1:num_obs
    
    rng(rng_fixed(i));

    seeds = rand(num_sim,num_v); % drawing random number from U[0,1]
    next_err = 99999*ones(size(seeds));
      
    next_err(:,1) = interp1(f1,ee1,seeds(:,1));
    next_err(:,2) = interp1(f2,ee2,seeds(:,2));   
    next_err(:,3) = interp1(f3,ee3,seeds(:,3));
    next_err(:,4) = interp1(f4,ee4,seeds(:,4));
    next_err(:,5) = interp1(f5,ee5,seeds(:,5));
    next_err(:,6) = interp1(f6,ee6,seeds(:,6));
    next_err(:,7) = interp1(f7,ee7,seeds(:,7));
    next_err(:,8) = interp1(f8,ee8,seeds(:,8));
    
    logasset_next = Elogasset_next(i,:) + next_err(:,1);
    equity_a_next = Eequity_a_next(i,:) + next_err(:,2);
    npf_a_next = Enpf_a_next(i,:) + next_err(:,3);
    realest_a_next = Erealest_a_next(i,:) + next_err(:,4);
    if aftersimulationfix == 1
        Jnpf = logical(npf_a_next <0);
        npf_a_next(Jnpf,:) = 0;
        Jrealest = logical(realest_a_next < 0);
        realest_a_next(Jrealest,:) = 0;
    end
	roa_next = Eroa_next(i,:) + next_err(:,5);
	grow_a_next = Egrow_a_next(i,:) + next_err(:,6);
    asset_next = Easset_next(i,:) + next_err(:,7);
    % asset_next = asset_next/(10^5); % NOT SURE if this is necessary
    unemp_state_next = Eunemp_next(i) + next_err(:,8);
    trend = yearq(i) - 102; %so that 1 represent yq(1985,4)
    trend_next = trend + 1;
    
        fspace = fundefn('spli',num_basis, min_equity_a,max_equity_a);
        equity_a_basis = funbas(fspace,equity_a_next); % note that we use "next" one here unlike logit code.
        equity_a_basis = full(equity_a_basis);
        equity_a_basis = equity_a_basis(:,1:end-1);
        % basis functions for npf_a;
        fspace = fundefn('spli',num_basis, min_npf_a,max_npf_a);
        npf_a_basis = funbas(fspace,npf_a_next);
        npf_a_basis = full(npf_a_basis);
        npf_a_basis = npf_a_basis(:,1:end-1);
        % basis functions for logasset;
        fspace = fundefn('spli',num_basis, min_logasset,max_logasset);
        logasset_basis = funbas(fspace,logasset_next);
        logasset_basis = full(logasset_basis);
        logasset_basis = logasset_basis(:,1:end-1);
        % basis functions for realEst_own_a;
        fspace = fundefn('spli',num_basis, min_realEst_own_a,max_realEst_own_a);
        realEst_own_a_basis = funbas(fspace,realest_a_next);
        realEst_own_a_basis = full(realEst_own_a_basis);
        realEst_own_a_basis = realEst_own_a_basis(:,1:end-1);
        % basis functions for roa_annual;
        fspace = fundefn('spli',num_basis, min_roa_annual,max_roa_annual);
        roa_annual_basis = funbas(fspace,roa_next);
        roa_annual_basis = full(roa_annual_basis);
        roa_annual_basis = roa_annual_basis(:,1:end-1);
        % basis functions for assGrowth_1y;
        fspace = fundefn('spli',num_basis, min_assGrowth_1y,max_assGrowth_1y);
        assGrowth_1y_basis = funbas(fspace,grow_a_next);
        assGrowth_1y_basis = full(assGrowth_1y_basis);
        assGrowth_1y_basis = assGrowth_1y_basis(:,1:end-1);
     %{
   [equity_a_all_basis,... %12(3*4) columns: [current spline (3columns), lag1 spline (3columns), lag2 ...]
		npf_a_all_basis,... %24
		logasset_all_basis,... %36
		realEst_own_a_all_basis,...  %48.
		roa_annual_all_basis,... %60
		assGrowth_1y_all_basis,... %72
        unemp_all, House, Senate, yearq, Bclose]; %73~76, 77, 78, 79, 80
   %}
    plogit_X = [equity_a_basis, repmat(equity_a_all_basis_temp(i,:),num_sim,1), ...
                npf_a_basis, repmat(npf_a_all_basis_temp(i,:),num_sim,1), ...
                logasset_basis, repmat(logasset_all_basis_temp(i,:),num_sim,1),...
                realEst_own_a_basis, repmat(realEst_own_a_all_basis_temp(i,:),num_sim,1),...
                roa_annual_basis, repmat(roa_annual_all_basis_temp(i,:),num_sim,1), ...
                assGrowth_1y_basis, repmat(assGrowth_1y_all_basis_temp(i,:),num_sim,1) ...
                unemp_state_next, repmat(unemp_all(i,1:end-1),num_sim,1),...
                repmat(House_next(i,:),num_sim,1), repmat(Senate_next(i,:),num_sim,1),ones(num_sim,1)];
    
    plogit_trend_X = [plogit_X(:,1:end-1), repmat(trend_next,num_sim,1), ones(num_sim,1)];
    
    prob_Bclose = exp(plogit_X*plogit_beta)./(1+exp(plogit_X*plogit_beta));
	Eprob_Bclose(i,:) = sum(prob_Bclose)/num_sim;
	ln_prob_Bclose = log(prob_Bclose);
	Eln_prob_Bclose(i,:) = sum(ln_prob_Bclose)/num_sim;
    %{
    %%% If the model with trend is used:
	prob_trend_Bclose = exp(plogit_trend_X*plogit_trend_beta)./(1+exp(plogit_trend_X*plogit_trend_beta));
	Eprob_trend_Bclose(i,:) = sum(prob_trend_Bclose)/num_sim;
	ln_prob_trend_Bclose = log(prob_trend_Bclose);
	Eln_prob_trend_Bclose(i,:) = sum(ln_prob_trend_Bclose)/num_sim;
    %}
    Xmc_next = [logasset_next, equity_a_next, npf_a_next, realest_a_next, roa_next, unemp_state_next, ones(num_sim,1)];
    mc_next_temp = Xmc_next*beta_mc;
    mc_next_temp = normcdf(mc_next_temp/sigma_mc,0,1).*mc_next_temp+sigma_mc*normpdf(mc_next_temp/sigma_mc,0,1);
    mc_next = mc_next_temp.*asset_next;
    Emc_next(i,:) = sum(mc_next)/num_sim; % already normalized values. (by /(10^5))
    J = logical(mc_next > 0);
	
    Enpf_a_next2(i,:) = sum(npf_a_next.*npf_a_next)/num_sim;
    Eroa_next2(i,:) = sum(roa_next.*roa_next)/num_sim;
    Elogasset_next2(i,:) = sum(logasset_next.*logasset_next)/num_sim;
end
toc
% takes about 450 seconds

IDyearqP = [BankID, yearq, Eprob_Bclose, Eln_prob_Bclose, Emc_next,...
    Enpf_a_next, Eroa_next, Enpf_a_next2, Eroa_next2, Egrow_a_next, Erealest_a_next, Elogasset_next, Elogasset_next2, ...
	Eequity_a_next, Eunemp_next];
%{
   %%% If the model with trend is used:
IDyearqP = [BankID, yearq, Eprob_Bclose, Eln_prob_Bclose, Emc_next,...
    Enpf_a_next, Eroa_next, Enpf_a_next2, Eroa_next2, Egrow_a_next, Erealest_a_next, Elogasset_next, Elogasset_next2, ...
	Eequity_a_next, Eunemp_next, Eprob_trend_Bclose, Eln_prob_trend_Bclose];
%}
save('C:\Users\domnisc1\Dropbox\structural econometrics course not shared\part 2\Banking_RFS15\output_files\ccpnext_param_sim5000_KLW.txt', 'IDyearqP', '-ASCII','-double','-tabs');

toc