% Variance of the estimates.

clear
load 'C:\Users\LoweryR\Documents\My Dropbox\UND_research_Banking\Final_files\output_files\dataforparam_KLW.mat';
load 'C:\Users\LoweryR\Documents\My Dropbox\UND_research_Banking\Final_files\output_files\theta_KLW.mat';

THETA = theta_KLW;
eulerconst = 0.577215665;
ONES = ones(size(Eprob_now,1),1);
Nobs = size(ONES,1);
Xpolinow = [House,Senate];
Xpolinext = [House_next,Senate_next];

% Problem with reporting Inf.
% The problem arises because one estimated probability is too close to 1.
% We assign it to be a bit bigger than second biggest one.

[aaa,bbb] = sortrows(Eprob_now); %aaa:sorted result, bbb:index
index_biggest = bbb(end);
index_2biggest = bbb(end-1);
Eprob_now(index_biggest,:) = Eprob_now(index_2biggest,:) + (1 - Eprob_now(index_2biggest,:))*0.9; 

beta = THETA(1);
sigma = THETA(2);
beta_NMC = THETA(3:10);
beta_NMC = beta_NMC';
Xnext = [ONES,Elogasset_next, Elogasset_next2, Enpf_a_next, Eroa_next, Erealest_a_next, Xpolinext];
Xnow = [ONES, logasset_all(:,1),logasset_all(:,1).*logasset_all(:,1), npf_a_all(:,1), roa_annual_all(:,1), realEst_own_a_all(:,1), Xpolinow];

Z = [ONES,logasset_all(:,1),logasset_all(:,1).*logasset_all(:,1),npf_a_all(:,1),roa_annual_all(:,1),Xpolinow(:,1),Xpolinow(:,2),equity_a_all(:,1),...
    realEst_own_a_all(:,1),assGrowth_1y_all(:,1),unemp_all(:,1)];
g1 = sigma*(log(1-Eprob_now)-log(Eprob_now)) + ...
    beta*sigma*Eln_prob_next - beta*(-Emc_next + Xnext*beta_NMC) + (-Emc_now + Xnow*beta_NMC);
g2 = g1.*Z(:,2);
g3 = g1.*Z(:,3);
g4 = g1.*Z(:,4);
g5 = g1.*Z(:,5);
g6 = g1.*Z(:,6);
g7 = g1.*Z(:,7); 
g8 = g1.*Z(:,8); %variable used in cost function
g9 = g1.*Z(:,9); 
g10 = g1.*Z(:,10); %variable used in transitions
g11 = g1.*Z(:,11);

m(1) = (1/Nobs)*sum(g1);
m(2) = (1/Nobs)*sum(g2);
m(3) = (1/Nobs)*sum(g3);
m(4) = (1/Nobs)*sum(g4);
m(5) = (1/Nobs)*sum(g5);
m(6) = (1/Nobs)*sum(g6);
m(7) = (1/Nobs)*sum(g7);
m(8) = (1/Nobs)*sum(g8);
m(9) = (1/Nobs)*sum(g9);
m(10) = (1/Nobs)*sum(g10);
m(11) = (1/Nobs)*sum(g11);

% Calculate AVar ignoring simulation and pre-estimation
% Hayashi p481: Avar = (G'WG)^-1 G'WSWG (G'WG)^-1
% Since it is exactly identified, Weighting matrix == I.
% So Avar = inv(G'G)G'SG inv(G'G)^-1
% G = E(dg/dtheta), S = E(g g') (Hayashi p481)
% But since G = [9 x 9]; Avar becomes inv(G)S inv(G')

g_all = [g1,g2,g3,g4,g5,g6,g7,g8,g9,g10,g11]';
Vhh = (1/Nobs)*(g_all*g_all');
% G = [N_g x N_theta]
% g1 = sigma*(log(1-Eprob_now)-log(Eprob_now)) + ...
%     beta*sigma*Eln_prob_next - beta*(-Emc_next + Xnext*beta_NMC) + (-Emc_now + Xnow*beta_NMC);

Hlambda = zeros(11,10);
Hlambda(1,1) = (1/Nobs)*sum(sigma*Eln_prob_next+Emc_next-Xnext*beta_NMC); %dg1/dbeta;
Hlambda(1,2) = (1/Nobs)*sum(log(1-Eprob_now)-log(Eprob_now)+beta*Eln_prob_next); %dg1/dsigma;
Hlambda(1,3) = (1/Nobs)*sum(-beta*Xnext(:,1)+Xnow(:,1));
Hlambda(1,4) = (1/Nobs)*sum(-beta*Xnext(:,2)+Xnow(:,2));
Hlambda(1,5) = (1/Nobs)*sum(-beta*Xnext(:,3)+Xnow(:,3));
Hlambda(1,6) = (1/Nobs)*sum(-beta*Xnext(:,4)+Xnow(:,4));
Hlambda(1,7) = (1/Nobs)*sum(-beta*Xnext(:,5)+Xnow(:,5));
Hlambda(1,8) = (1/Nobs)*sum(-beta*Xnext(:,6)+Xnow(:,6));
Hlambda(1,9) = (1/Nobs)*sum(-beta*Xnext(:,7)+Xnow(:,7));
Hlambda(1,10) = (1/Nobs)*sum(-beta*Xnext(:,8)+Xnow(:,8));

Hlambda(2,1) = (1/Nobs)*sum((sigma*Eln_prob_next+Emc_next-Xnext*beta_NMC).*Z(:,2)); %dg1/dbeta;
Hlambda(2,2) = (1/Nobs)*sum((log(1-Eprob_now)-log(Eprob_now)+beta*Eln_prob_next).*Z(:,2)); %dg1/dsigma;
Hlambda(2,3) = (1/Nobs)*sum((-beta*Xnext(:,1)+Xnow(:,1)).*Z(:,2));
Hlambda(2,4) = (1/Nobs)*sum((-beta*Xnext(:,2)+Xnow(:,2)).*Z(:,2));
Hlambda(2,5) = (1/Nobs)*sum((-beta*Xnext(:,3)+Xnow(:,3)).*Z(:,2));
Hlambda(2,6) = (1/Nobs)*sum((-beta*Xnext(:,4)+Xnow(:,4)).*Z(:,2));
Hlambda(2,7) = (1/Nobs)*sum((-beta*Xnext(:,5)+Xnow(:,5)).*Z(:,2));
Hlambda(2,8) = (1/Nobs)*sum((-beta*Xnext(:,6)+Xnow(:,6)).*Z(:,2));
Hlambda(2,9) = (1/Nobs)*sum((-beta*Xnext(:,7)+Xnow(:,7)).*Z(:,2));
Hlambda(2,10) = (1/Nobs)*sum((-beta*Xnext(:,8)+Xnow(:,8)).*Z(:,2));

Hlambda(3,1) = (1/Nobs)*sum((sigma*Eln_prob_next+Emc_next-Xnext*beta_NMC).*Z(:,3)); %dg1/dbeta;
Hlambda(3,2) = (1/Nobs)*sum((log(1-Eprob_now)-log(Eprob_now)+beta*Eln_prob_next).*Z(:,3)); %dg1/dsigma;
Hlambda(3,3) = (1/Nobs)*sum((-beta*Xnext(:,1)+Xnow(:,1)).*Z(:,3));
Hlambda(3,4) = (1/Nobs)*sum((-beta*Xnext(:,2)+Xnow(:,2)).*Z(:,3));
Hlambda(3,5) = (1/Nobs)*sum((-beta*Xnext(:,3)+Xnow(:,3)).*Z(:,3));
Hlambda(3,6) = (1/Nobs)*sum((-beta*Xnext(:,4)+Xnow(:,4)).*Z(:,3));
Hlambda(3,7) = (1/Nobs)*sum((-beta*Xnext(:,5)+Xnow(:,5)).*Z(:,3));
Hlambda(3,8) = (1/Nobs)*sum((-beta*Xnext(:,6)+Xnow(:,6)).*Z(:,3));
Hlambda(3,9) = (1/Nobs)*sum((-beta*Xnext(:,7)+Xnow(:,7)).*Z(:,3));
Hlambda(3,10) = (1/Nobs)*sum((-beta*Xnext(:,8)+Xnow(:,8)).*Z(:,3));

Hlambda(4,1) = (1/Nobs)*sum((sigma*Eln_prob_next+Emc_next-Xnext*beta_NMC).*Z(:,4)); %dg1/dbeta;
Hlambda(4,2) = (1/Nobs)*sum((log(1-Eprob_now)-log(Eprob_now)+beta*Eln_prob_next).*Z(:,4)); %dg1/dsigma;
Hlambda(4,3) = (1/Nobs)*sum((-beta*Xnext(:,1)+Xnow(:,1)).*Z(:,4));
Hlambda(4,4) = (1/Nobs)*sum((-beta*Xnext(:,2)+Xnow(:,2)).*Z(:,4));
Hlambda(4,5) = (1/Nobs)*sum((-beta*Xnext(:,3)+Xnow(:,3)).*Z(:,4));
Hlambda(4,6) = (1/Nobs)*sum((-beta*Xnext(:,4)+Xnow(:,4)).*Z(:,4));
Hlambda(4,7) = (1/Nobs)*sum((-beta*Xnext(:,5)+Xnow(:,5)).*Z(:,4));
Hlambda(4,8) = (1/Nobs)*sum((-beta*Xnext(:,6)+Xnow(:,6)).*Z(:,4));
Hlambda(4,9) = (1/Nobs)*sum((-beta*Xnext(:,7)+Xnow(:,7)).*Z(:,4));
Hlambda(4,10) = (1/Nobs)*sum((-beta*Xnext(:,8)+Xnow(:,8)).*Z(:,4));

Hlambda(5,1) = (1/Nobs)*sum((sigma*Eln_prob_next+Emc_next-Xnext*beta_NMC).*Z(:,5)); %dg1/dbeta;
Hlambda(5,2) = (1/Nobs)*sum((log(1-Eprob_now)-log(Eprob_now)+beta*Eln_prob_next).*Z(:,5)); %dg1/dsigma;
Hlambda(5,3) = (1/Nobs)*sum((-beta*Xnext(:,1)+Xnow(:,1)).*Z(:,5));
Hlambda(5,4) = (1/Nobs)*sum((-beta*Xnext(:,2)+Xnow(:,2)).*Z(:,5));
Hlambda(5,5) = (1/Nobs)*sum((-beta*Xnext(:,3)+Xnow(:,3)).*Z(:,5));
Hlambda(5,6) = (1/Nobs)*sum((-beta*Xnext(:,4)+Xnow(:,4)).*Z(:,5));
Hlambda(5,7) = (1/Nobs)*sum((-beta*Xnext(:,5)+Xnow(:,5)).*Z(:,5));
Hlambda(5,8) = (1/Nobs)*sum((-beta*Xnext(:,6)+Xnow(:,6)).*Z(:,5));
Hlambda(5,9) = (1/Nobs)*sum((-beta*Xnext(:,7)+Xnow(:,7)).*Z(:,5));
Hlambda(5,10) = (1/Nobs)*sum((-beta*Xnext(:,8)+Xnow(:,8)).*Z(:,5));

Hlambda(6,1) = (1/Nobs)*sum((sigma*Eln_prob_next+Emc_next-Xnext*beta_NMC).*Z(:,6)); %dg1/dbeta;
Hlambda(6,2) = (1/Nobs)*sum((log(1-Eprob_now)-log(Eprob_now)+beta*Eln_prob_next).*Z(:,6)); %dg1/dsigma;
Hlambda(6,3) = (1/Nobs)*sum((-beta*Xnext(:,1)+Xnow(:,1)).*Z(:,6));
Hlambda(6,4) = (1/Nobs)*sum((-beta*Xnext(:,2)+Xnow(:,2)).*Z(:,6));
Hlambda(6,5) = (1/Nobs)*sum((-beta*Xnext(:,3)+Xnow(:,3)).*Z(:,6));
Hlambda(6,6) = (1/Nobs)*sum((-beta*Xnext(:,4)+Xnow(:,4)).*Z(:,6));
Hlambda(6,7) = (1/Nobs)*sum((-beta*Xnext(:,5)+Xnow(:,5)).*Z(:,6));
Hlambda(6,8) = (1/Nobs)*sum((-beta*Xnext(:,6)+Xnow(:,6)).*Z(:,6));
Hlambda(6,9) = (1/Nobs)*sum((-beta*Xnext(:,7)+Xnow(:,7)).*Z(:,6));
Hlambda(6,10) = (1/Nobs)*sum((-beta*Xnext(:,8)+Xnow(:,8)).*Z(:,6));

Hlambda(7,1) = (1/Nobs)*sum((sigma*Eln_prob_next+Emc_next-Xnext*beta_NMC).*Z(:,7)); %dg1/dbeta;
Hlambda(7,2) = (1/Nobs)*sum((log(1-Eprob_now)-log(Eprob_now)+beta*Eln_prob_next).*Z(:,7)); %dg1/dsigma;
Hlambda(7,3) = (1/Nobs)*sum((-beta*Xnext(:,1)+Xnow(:,1)).*Z(:,7));
Hlambda(7,4) = (1/Nobs)*sum((-beta*Xnext(:,2)+Xnow(:,2)).*Z(:,7));
Hlambda(7,5) = (1/Nobs)*sum((-beta*Xnext(:,3)+Xnow(:,3)).*Z(:,7));
Hlambda(7,6) = (1/Nobs)*sum((-beta*Xnext(:,4)+Xnow(:,4)).*Z(:,7));
Hlambda(7,7) = (1/Nobs)*sum((-beta*Xnext(:,5)+Xnow(:,5)).*Z(:,7));
Hlambda(7,8) = (1/Nobs)*sum((-beta*Xnext(:,6)+Xnow(:,6)).*Z(:,7));
Hlambda(7,9) = (1/Nobs)*sum((-beta*Xnext(:,7)+Xnow(:,7)).*Z(:,7));
Hlambda(7,10) = (1/Nobs)*sum((-beta*Xnext(:,8)+Xnow(:,8)).*Z(:,7));

Hlambda(8,1) = (1/Nobs)*sum((sigma*Eln_prob_next+Emc_next-Xnext*beta_NMC).*Z(:,8)); %dg1/dbeta;
Hlambda(8,2) = (1/Nobs)*sum((log(1-Eprob_now)-log(Eprob_now)+beta*Eln_prob_next).*Z(:,8)); %dg1/dsigma;
Hlambda(8,3) = (1/Nobs)*sum((-beta*Xnext(:,1)+Xnow(:,1)).*Z(:,8));
Hlambda(8,4) = (1/Nobs)*sum((-beta*Xnext(:,2)+Xnow(:,2)).*Z(:,8));
Hlambda(8,5) = (1/Nobs)*sum((-beta*Xnext(:,3)+Xnow(:,3)).*Z(:,8));
Hlambda(8,6) = (1/Nobs)*sum((-beta*Xnext(:,4)+Xnow(:,4)).*Z(:,8));
Hlambda(8,7) = (1/Nobs)*sum((-beta*Xnext(:,5)+Xnow(:,5)).*Z(:,8));
Hlambda(8,8) = (1/Nobs)*sum((-beta*Xnext(:,6)+Xnow(:,6)).*Z(:,8));
Hlambda(8,9) = (1/Nobs)*sum((-beta*Xnext(:,7)+Xnow(:,7)).*Z(:,8));
Hlambda(8,10) = (1/Nobs)*sum((-beta*Xnext(:,8)+Xnow(:,8)).*Z(:,8));

Hlambda(9,1) = (1/Nobs)*sum((sigma*Eln_prob_next+Emc_next-Xnext*beta_NMC).*Z(:,9)); %dg1/dbeta;
Hlambda(9,2) = (1/Nobs)*sum((log(1-Eprob_now)-log(Eprob_now)+beta*Eln_prob_next).*Z(:,9)); %dg1/dsigma;
Hlambda(9,3) = (1/Nobs)*sum((-beta*Xnext(:,1)+Xnow(:,1)).*Z(:,9));
Hlambda(9,4) = (1/Nobs)*sum((-beta*Xnext(:,2)+Xnow(:,2)).*Z(:,9));
Hlambda(9,5) = (1/Nobs)*sum((-beta*Xnext(:,3)+Xnow(:,3)).*Z(:,9));
Hlambda(9,6) = (1/Nobs)*sum((-beta*Xnext(:,4)+Xnow(:,4)).*Z(:,9));
Hlambda(9,7) = (1/Nobs)*sum((-beta*Xnext(:,5)+Xnow(:,5)).*Z(:,9));
Hlambda(9,8) = (1/Nobs)*sum((-beta*Xnext(:,6)+Xnow(:,6)).*Z(:,9));
Hlambda(9,9) = (1/Nobs)*sum((-beta*Xnext(:,7)+Xnow(:,7)).*Z(:,9));
Hlambda(9,10) = (1/Nobs)*sum((-beta*Xnext(:,8)+Xnow(:,8)).*Z(:,9));

Hlambda(10,1) = (1/Nobs)*sum((sigma*Eln_prob_next+Emc_next-Xnext*beta_NMC).*Z(:,10)); %dg1/dbeta;
Hlambda(10,2) = (1/Nobs)*sum((log(1-Eprob_now)-log(Eprob_now)+beta*Eln_prob_next).*Z(:,10)); %dg1/dsigma;
Hlambda(10,3) = (1/Nobs)*sum((-beta*Xnext(:,1)+Xnow(:,1)).*Z(:,10));
Hlambda(10,4) = (1/Nobs)*sum((-beta*Xnext(:,2)+Xnow(:,2)).*Z(:,10));
Hlambda(10,5) = (1/Nobs)*sum((-beta*Xnext(:,3)+Xnow(:,3)).*Z(:,10));
Hlambda(10,6) = (1/Nobs)*sum((-beta*Xnext(:,4)+Xnow(:,4)).*Z(:,10));
Hlambda(10,7) = (1/Nobs)*sum((-beta*Xnext(:,5)+Xnow(:,5)).*Z(:,10));
Hlambda(10,8) = (1/Nobs)*sum((-beta*Xnext(:,6)+Xnow(:,6)).*Z(:,10));
Hlambda(10,9) = (1/Nobs)*sum((-beta*Xnext(:,7)+Xnow(:,7)).*Z(:,10));
Hlambda(10,10) = (1/Nobs)*sum((-beta*Xnext(:,8)+Xnow(:,8)).*Z(:,10));

Hlambda(11,1) = (1/Nobs)*sum((sigma*Eln_prob_next+Emc_next-Xnext*beta_NMC).*Z(:,11)); %dg1/dbeta;
Hlambda(11,2) = (1/Nobs)*sum((log(1-Eprob_now)-log(Eprob_now)+beta*Eln_prob_next).*Z(:,11)); %dg1/dsigma;
Hlambda(11,3) = (1/Nobs)*sum((-beta*Xnext(:,1)+Xnow(:,1)).*Z(:,11));
Hlambda(11,4) = (1/Nobs)*sum((-beta*Xnext(:,2)+Xnow(:,2)).*Z(:,11));
Hlambda(11,5) = (1/Nobs)*sum((-beta*Xnext(:,3)+Xnow(:,3)).*Z(:,11));
Hlambda(11,6) = (1/Nobs)*sum((-beta*Xnext(:,4)+Xnow(:,4)).*Z(:,11));
Hlambda(11,7) = (1/Nobs)*sum((-beta*Xnext(:,5)+Xnow(:,5)).*Z(:,11));
Hlambda(11,8) = (1/Nobs)*sum((-beta*Xnext(:,6)+Xnow(:,6)).*Z(:,11));
Hlambda(11,9) = (1/Nobs)*sum((-beta*Xnext(:,7)+Xnow(:,7)).*Z(:,11));
Hlambda(11,10) = (1/Nobs)*sum((-beta*Xnext(:,8)+Xnow(:,8)).*Z(:,11));
Hlambda(11,10) = (1/Nobs)*sum((-beta*Xnext(:,8)+Xnow(:,8)).*Z(:,11));

% Avar = inv(Hlambda)*Vhh*inv(Hlambda'); % WHEN #moments == #param
A = Vhh\Hlambda;
Avar = inv(Hlambda'*A); % WHEN #moments > #param (Hayashi p481)
THETA_var =  (1/Nobs)*Avar;
SD_var = sqrt(diag(THETA_var));
T_noadjust = [THETA',SD_var,THETA'./SD_var];

Avar_lambda = Avar;

% Let's take care of pre-estimation.
% Adjust pre-estimation of cost function
load 'C:\Users\LoweryR\Documents\My Dropbox\UND_research_Banking\Final_files\output_files\ccpnext_param_sim5000_KLW.txt';
%{
'ccpnext_param_sim500_KLW.mat' contains ...
    'Enpf_next','Eequity_next','Eequity2_l_next','Eunemp_state_liab_next',...
	'Edenomi','Eequity_l_basis_denomi','Enpf_l_basis_denomi','Enet_inc_l_basis_denomi','Eliability_basis_denomi');
%}
beta_mc_sigma = xlsread('C:\Users\LoweryR\Documents\My Dropbox\UND_research_Banking\Final_files\data_files\coeff_resid_KLW.xlsx','cost_beta_v26');
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(Nobs,1)];
Xmc_next = [Elogasset_next,Eequity_a_next,Enpf_a_next,Erealest_a_next,Eroa_next,Eunemp_next, ones(Nobs,1)];
mc_now_temp = Xmc_now*beta_mc;
mc_next_temp = Xmc_next*beta_mc;
nvar_cost = length(Xmc_now(1,:));
% Note: Cost function is Tobit.  So dE(y|x)dparam = normcdf*x See Wooldridge p523 Eq.16.16.(2002 edition)
% g1 = sigma*(log(1-Eprob_now)-log(Eprob_now)) + ...
%     beta*sigma*Eln_prob_next - beta*(-Emc_next + Xnext*beta_NMC) + (-Emc_now + Xnow*beta_NMC);
dg1dcost = beta*Xmc_next.*repmat(normcdf(mc_next_temp/sigma_mc,0,1).*Easset_next,1,nvar_cost) ...
    - Xmc_now.*repmat(normcdf(mc_now_temp/sigma_mc,0,1).*asset_all(:,1),1,nvar_cost); %[Nobs x 5]
dg2dcost = dg1dcost.*repmat(Z(:,2),1,nvar_cost);
dg3dcost = dg1dcost.*repmat(Z(:,3),1,nvar_cost);
dg4dcost = dg1dcost.*repmat(Z(:,4),1,nvar_cost);
dg5dcost = dg1dcost.*repmat(Z(:,5),1,nvar_cost);
dg6dcost = dg1dcost.*repmat(Z(:,6),1,nvar_cost);
dg7dcost = dg1dcost.*repmat(Z(:,7),1,nvar_cost);
dg8dcost = dg1dcost.*repmat(Z(:,8),1,nvar_cost);
dg9dcost = dg1dcost.*repmat(Z(:,9),1,nvar_cost);
dg10dcost = dg1dcost.*repmat(Z(:,10),1,nvar_cost);
dg11dcost = dg1dcost.*repmat(Z(:,11),1,nvar_cost);

Hcost = zeros(11,nvar_cost);
Hcost(1,:) = (1/Nobs)*sum(dg1dcost);
Hcost(2,:) = (1/Nobs)*sum(dg2dcost);
Hcost(3,:) = (1/Nobs)*sum(dg3dcost);
Hcost(4,:) = (1/Nobs)*sum(dg4dcost);
Hcost(5,:) = (1/Nobs)*sum(dg5dcost);
Hcost(6,:) = (1/Nobs)*sum(dg6dcost);
Hcost(7,:) = (1/Nobs)*sum(dg7dcost);
Hcost(8,:) = (1/Nobs)*sum(dg8dcost);
Hcost(9,:) = (1/Nobs)*sum(dg9dcost);
Hcost(10,:) = (1/Nobs)*sum(dg10dcost);
Hcost(11,:) = (1/Nobs)*sum(dg11dcost);

% Let's take care of pre-estimation.
% Adjust pre-estimation of logit function
load 'C:\Users\LoweryR\Documents\My Dropbox\UND_research_Banking\Final_files\data_files\dataforlogit_KLW.txt';
trend = dataforlogit_KLW(:,end-1) - 102; %v79 is the trend, v80(the last one) is the Bclose
Nobs_logit = size(trend,1);
Xlogit_now = [dataforlogit_KLW(:,1:end-2), ones(Nobs_logit,1)];
%{
dataforlogit_KLW.txt has 80 columns.
[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
%}
num_lags = 4; % number of sequences in each variable; here the last 4th lag is not included so curent, lag1,..lag3 => 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 colums to avoid multicolinearity (remember: sum(basis,2) == 1)
num_l_b = num_lags*num_basis_a;
equity_part = dataforlogit_KLW(:,1:num_l_b);
npf_part = dataforlogit_KLW(:,1+num_l_b:2*num_l_b);
logasset_part = dataforlogit_KLW(:,2*num_l_b+1:3*num_l_b);
realEst_part = dataforlogit_KLW(:,3*num_l_b+1:4*num_l_b);
roa_part = dataforlogit_KLW(:,4*num_l_b+1:5*num_l_b);
assGrowth_part = dataforlogit_KLW(:,5*num_l_b+1:6*num_l_b);
unemp_part = dataforlogit_KLW(:,6*num_l_b+1:6*num_l_b+4);

% NOTE! we use min/max out of ALL samples. The knot locations are consistent accorss 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));

% basis functions for future equity_a;
fspace = fundefn('spli',num_basis, min_equity_a, max_equity_a);
equity_a_basis = funbas(fspace,Eequity_a_next);
equity_a_basis = full(equity_a_basis);
equity_next_basis = equity_a_basis(:,1:end-1);
% basis functions for future npf_a;
fspace = fundefn('spli',num_basis, min_npf_a,max_npf_a);
npf_a_basis = funbas(fspace,Enpf_a_next);
npf_a_basis = full(npf_a_basis);
npf_next_basis = npf_a_basis(:,1:end-1);
% basis functions for future logasset;
fspace = fundefn('spli',num_basis, min_logasset,max_logasset);
logasset_basis = funbas(fspace,Elogasset_next);
logasset_basis = full(logasset_basis);
logasset_next_basis = logasset_basis(:,1:end-1);
% basis functions for future realEst_own_a;
fspace = fundefn('spli',num_basis, min_realEst_own_a,max_realEst_own_a);
realEst_own_a_basis = funbas(fspace,Erealest_a_next);
realEst_own_a_basis = full(realEst_own_a_basis);
realEst_next_basis = realEst_own_a_basis(:,1:end-1);
% basis functions for future roa_annual;
fspace = fundefn('spli',num_basis, min_roa_annual,max_roa_annual);
roa_annual_basis = funbas(fspace,Eroa_next);
roa_annual_basis = full(roa_annual_basis);
roa_next_basis = roa_annual_basis(:,1:end-1);
% basis functions for future assGrowth_1y;
fspace = fundefn('spli',num_basis, min_assGrowth_1y,max_assGrowth_1y);
assGrowth_1y_basis = funbas(fspace,Egrow_a_next);
assGrowth_1y_basis = full(assGrowth_1y_basis);
assGrowth_next_basis = assGrowth_1y_basis(:,1:end-1);

Xlogit_next = [equity_next_basis, equity_part(:,1:(num_lags-1)*num_basis_a),... % consist of future basis, and the basis from the one period forwarded basis
    npf_next_basis, npf_part(:,1:(num_lags-1)*num_basis_a),...
    logasset_next_basis, logasset_part(:,1:(num_lags-1)*num_basis_a),...
    realEst_next_basis, realEst_part(:,1:(num_lags-1)*num_basis_a),...
    roa_next_basis, roa_part(:,1:(num_lags-1)*num_basis_a),...
    assGrowth_next_basis, assGrowth_part(:,1:(num_lags-1)*num_basis_a),...
    Eunemp_next, unemp_part(:,1:end-1),...
    House_next, Senate_next,ones(Nobs_logit,1)];
nvar_logit = size(Xlogit_now,2);
dg1dlogit = -THETA(2)*Xlogit_now + THETA(1)*THETA(2)*repmat(1-Eprob_next,1,nvar_logit).*Xlogit_next; %[Nobs x 5]
dg2dlogit = dg1dlogit.*repmat(Z(:,2),1,nvar_logit);
dg3dlogit = dg1dlogit.*repmat(Z(:,3),1,nvar_logit);
dg4dlogit = dg1dlogit.*repmat(Z(:,4),1,nvar_logit);
dg5dlogit = dg1dlogit.*repmat(Z(:,5),1,nvar_logit);
dg6dlogit = dg1dlogit.*repmat(Z(:,6),1,nvar_logit);
dg7dlogit = dg1dlogit.*repmat(Z(:,7),1,nvar_logit);
dg8dlogit = dg1dlogit.*repmat(Z(:,8),1,nvar_logit);
dg9dlogit = dg1dlogit.*repmat(Z(:,9),1,nvar_logit);
dg10dlogit = dg1dlogit.*repmat(Z(:,10),1,nvar_logit);
dg11dlogit = dg1dlogit.*repmat(Z(:,11),1,nvar_logit);

Hlogit = zeros(11,nvar_logit);
Hlogit(1,:) = (1/Nobs)*sum(dg1dlogit);
Hlogit(2,:) = (1/Nobs)*sum(dg2dlogit);
Hlogit(3,:) = (1/Nobs)*sum(dg3dlogit);
Hlogit(4,:) = (1/Nobs)*sum(dg4dlogit);
Hlogit(5,:) = (1/Nobs)*sum(dg5dlogit);
Hlogit(6,:) = (1/Nobs)*sum(dg6dlogit);
Hlogit(7,:) = (1/Nobs)*sum(dg7dlogit);
Hlogit(8,:) = (1/Nobs)*sum(dg8dlogit);
Hlogit(9,:) = (1/Nobs)*sum(dg9dlogit);
Hlogit(10,:) = (1/Nobs)*sum(dg10dlogit);
Hlogit(11,:) = (1/Nobs)*sum(dg11dlogit);

load 'C:\Users\LoweryR\Documents\My Dropbox\UND_research_Banking\Final_files\output_files\logit_covb_KLW_nonames.txt';% note that you have costANDtransition_v28v29.do ...
% Note that "_nonames" means the first row and first column of variable names are deleted so that Matlab can read the file.
Avar_logit = Nobs_logit*logit_covb_KLW_nonames;
load 'C:\Users\LoweryR\Documents\My Dropbox\UND_research_Banking\Final_files\data_files\cost_covb_KLW_nonames.txt';
% note that "_nonames" means the first row and first column of variable names are deleted so that Matlab
% can read the file.  Also for this file, the last column and the last
% numeric row are removed since those are correlation with "sigma".
Nobs_cost = 3678; % I think we should use the number of observations actually used to estimate cost fn.
Avar_cost = Nobs_cost*cost_covb_KLW_nonames;

temp_matrix = Vhh\Hlambda;
Glambda_W_Glambda = Hlambda'*temp_matrix;
temp_matrix = Vhh\Hcost;
Glambda_W_Gcost = Hlambda'*temp_matrix; 
temp_matrix = Glambda_W_Gcost*Avar_cost*Glambda_W_Gcost';
Avar_adj_bycost = Avar_lambda*temp_matrix*Avar_lambda;
Avar_adj_bycost_test = Avar_lambda*Glambda_W_Gcost*Avar_cost*Glambda_W_Gcost'*Avar_lambda;

temp_matrix = Glambda_W_Glambda\Glambda_W_Gcost;
Avar_adj_bycost_temp = temp_matrix*Avar_cost*Glambda_W_Gcost'*Avar_lambda;

temp_matrix = Vhh\Hlogit;
Glambda_W_Glogit = Hlambda'*temp_matrix; 
temp_matrix = Glambda_W_Glogit*Avar_logit*Glambda_W_Glogit';
Avar_adj_bylogit = Avar_lambda*temp_matrix*Avar_lambda;

Avar_adj_bylogit_test = Avar_lambda*Glambda_W_Glogit*Avar_logit*Glambda_W_Glogit'*Avar_lambda;

temp_matrix = Glambda_W_Glambda\Glambda_W_Glogit;
Avar_adj_bylogit_temp = temp_matrix*Avar_logit*Glambda_W_Glogit'*Avar_lambda;


% Avar = Avar_lambda + inv(Hlambda)*Hcost*(Avar_cost)*Hcost'*(inv(Hlambda))' + inv(Hlambda)*Hlogit*(Avar_logit)*Hlogit'*(inv(Hlambda))';
Avar = Avar_lambda + Avar_adj_bycost + Avar_adj_bylogit;
THETA_var =  (1/Nobs)*Avar;
SD_var = sqrt(diag(THETA_var));
T_adjust = [THETA',SD_var,THETA'./SD_var];

% This one uses "Nobs" for all pre-estimated Avar.
Avar_logit_A = Nobs*logit_covb_KLW_nonames;
Avar_cost_A = Nobs*cost_covb_KLW_nonames;

temp_matrix = Vhh\Hlambda;
Glambda_W_Glambda = Hlambda'*temp_matrix;
temp_matrix = Vhh\Hcost;
Glambda_W_Gcost = Hlambda'*temp_matrix; 
temp_matrix = Glambda_W_Gcost*Avar_cost_A*Glambda_W_Gcost';
Avar_adj_bycost_A = Avar_lambda*temp_matrix*Avar_lambda;
temp_matrix = Glambda_W_Glambda\Glambda_W_Gcost;
Avar_adj_bycost_A_temp = temp_matrix*Avar_cost_A*Glambda_W_Gcost'*Avar_lambda;

temp_matrix = Vhh\Hlogit;
Glambda_W_Glogit = Hlambda'*temp_matrix; 
temp_matrix = Glambda_W_Glogit*Avar_logit_A*Glambda_W_Glogit';
Avar_adj_bylogit_A = Avar_lambda*temp_matrix*Avar_lambda;
temp_matrix = Glambda_W_Glambda\Glambda_W_Glogit;
Avar_adj_bylogit_A_temp = temp_matrix*Avar_logit_A*Glambda_W_Glogit'*Avar_lambda;

% Avar_A = Avar_lambda + inv(Hlambda)*Hcost*(Avar_cost_A)*Hcost'*(inv(Hlambda))' + inv(Hlambda)*Hlogit*(Avar_logit_A)*Hlogit'*(inv(Hlambda))';
Avar_A = Avar_lambda + Avar_adj_bycost_A + Avar_adj_bylogit_A;
THETA_var_A =  (1/Nobs)*Avar_A;
SD_var_A = sqrt(diag(THETA_var_A));
T_adjust_A = [THETA',SD_var_A,THETA'./SD_var_A]; % parameters: beta, sigma, constant,... 

save('C:\Users\LoweryR\Documents\My Dropbox\UND_research_Banking\Final_files\output_files\T_KLW.mat','T_adjust_A'); % The results reported in Table 7, Model 1.