
function  [THETA_new,f_old,f_new,minf,THETA_minf] = fn_obj_KLW(THETA)

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

%{
'dataforparam_KLW.mat' variables:
    'BankID', 'yearq', 'Bclose',...
'Eprob_now','Eln_prob_next','Eprob_next','Eln_prob_t_next','Eprob_t_next','Emc_now','Emc_next',...
'equity_a_all', 'npf_a_all', 'logasset_all', 'realEst_own_a_all', 'roa_annual_all', 'unemp_all',...
'assGrowth_1y_all', 'House', 'Senate', 'House_next', 'Senate_next', 'asset_all',...
'Elogasset_next', 'Elogasset_next2','Eequity_a_next', 'Enpf_a_next', 'Enpf_a_next2', 'Erealest_a_next', 'Eroa_next', 'Eroa_next2', 'Egrow_a_next', 'Easset_next','Eunemp_next');
%}

% 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);

W = (1/Nobs)*(Z'*Z);
f_temp = W\m';
f_old = m*f_temp;

% See what happens if one of parameters moves around while others are fixed.

num_range = 100;
extend = 1;
portion = 0.01;
THETA_range = [THETA(1) - THETA(1)*portion*num_range:THETA(1)*portion:THETA(1) + extend*THETA(1)*portion*num_range];
f_range = NaN(length(THETA),length(THETA_range));
x_range = NaN(length(THETA),length(THETA_range));

for i_param = 1:length(THETA)
    THETA_range = [THETA(i_param) - THETA(i_param)*portion*num_range:THETA(i_param)*portion:THETA(i_param) + extend*THETA(i_param)*portion*num_range];
    for k = 1:length(THETA_range)
        if i_param == 1
            g1 =  THETA(2)*(log(1-Eprob_now)-log(Eprob_now)) + ...
                THETA_range(k)*THETA(2)*Eln_prob_next - THETA_range(k)*(-Emc_next + Xnext*THETA(3:10)') + (-Emc_now + Xnow*THETA(3:10)');
        elseif i_param == 2
            g1 = THETA_range(k)*(log(1-Eprob_now)-log(Eprob_now)) + ...
                THETA(1)*THETA_range(k)*Eln_prob_next - THETA(1)*(-Emc_next + Xnext*THETA(3:10)') + (-Emc_now + Xnow*THETA(3:10)');
        end
        for kk = 3:length(THETA)
            if i_param == kk
                beta_NMC_temp = THETA(3:10)';
                beta_NMC_temp(kk-2) = THETA_range(k);
                g1 = THETA(2)*(log(1-Eprob_now)-log(Eprob_now)) + ...
                    THETA(1)*THETA(2)*Eln_prob_next - THETA(1)*(-Emc_next + Xnext*beta_NMC_temp) + (-Emc_now + Xnow*beta_NMC_temp);
            end
        end
        
        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);
        
        W = (1/Nobs)*(Z'*Z);
        
        f_temp = W\m';
        f_range(i_param,k) = m*f_temp;
    end
    x_range(i_param,:) = THETA_range;
end

minf = min(min(f_range));
[theparam,theloc] = find(f_range == minf); 
THETA_minf = zeros(size(THETA));
if length(theparam) == 1
    THETA_minf = THETA;
    THETA_minf(theparam) = x_range(theparam,theloc);
end

THETA_new = NaN(size(THETA));
for i = 1:length(THETA)
    J = find(f_range(i,:) == min(f_range(i,:)));
    THETA_new(i) = x_range(i,J);
end

beta_new = THETA_new(1);
sigma_new = THETA_new(2);
beta_NMC_new = THETA_new(3:10);
beta_NMC_new = beta_NMC_new';
g1 = sigma_new*(log(1-Eprob_now)-log(Eprob_now)) + ...
    beta_new*sigma_new*Eln_prob_next - beta_new*(-Emc_next + Xnext*beta_NMC_new) + (-Emc_now + Xnow*beta_NMC_new);
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);

W = (1/Nobs)*(Z'*Z);

f_temp = W\m';
f_new = m*f_temp;