clear all;
cd 'C:\Documents and Settings\mw109\Desktop\teaching\applied econometrics\matlab_applied_micro';

format compact;
diary amicro_num_methods_diary

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Example 1: Numerical Derivatives
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

display('Numerical Derivatives Example');;

%approximate x.^2 using finite difference

h = 0.01;
X0 = 10;

%one sided
deriv_est1 = (1/(h * X0)) * ( (X0 + (h * X0) ).^3 - (X0).^3)

%two sided
deriv_est2 = (1/(2 * h * X0)) * ( (X0 + (h * X0) ).^3 - (X0 - (h * X0)).^3)

actual = 3 * X0^2


%using Miranda and Fackler h

%one sided
h = max(abs(X0),1) * sqrt(eps)
deriv_est1 = (1/(h * X0)) * ( (X0 + (h * X0) ).^3 - (X0).^3)

%two sided
h = max(abs(X0),1) * eps.^(1/3)
deriv_est2 = (1/(2 * h * X0)) * ( (X0 + (h * X0) ).^3 - (X0 - (h * X0)).^3)




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Example 2: Gauss-Hermite Quadrature
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;

for order = 1:2
EY(order) =  pi^(-1/2) * ( 0.886 * (sqrt(2) * 0.707)^order ...
+ 0.886 * (sqrt(2) * (-0.707))^order )
end;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Example 3: Newton-Raphson
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;

display('Newton-Raphson Example');;

%plot function
theta_vals = -10:0.1:10;
for j = 1:length(theta_vals)
    funct_evals(j) = newton_func(theta_vals(j));
end;

plot(theta_vals,funct_evals,theta_vals,zeros(length(theta_vals),1))


%estimation options for simplex routine    
options = optimset('LargeScale','off','TolFun', 0.0001, 'TolX', 0.0001, 'MaxIter',100000, ...
  'Display','iter','MaxFunEvals',5000); 


param_start = 3;

%solve using simplex
[param_estimates,fval,exitflag] = fminsearch('newton_func',param_start,options)

%solve using newton-raphson (with additional Matlab features)
[param_estimates,fval,exitflag] = fminunc('newton_func',param_start,options)


%solve using simple newton-raphson

param_old = param_start;
iter_max = 100;
crit = 100;
iter = 0;
while crit > 0.001;

iter = iter + 1    
if iter > iter_max;
    break;
else
    
G = -4 * param_old + param_old^2
H = -4 + 2 * param_old
param_new = param_old - G * inv(H)
crit = abs((param_new - param_old) / param_old)

param_old = param_new;

end; end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Example 4: Simplex
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;

display('Simplex Example');;

%starting parameters;
param_start = [1/2; 1/2]
J = 2;

%solve using Matlab's fminsearch
%estimation options for simplex routine    
options = optimset('LargeScale','off','TolFun', 0.0001, 'TolX', 0.0001, 'MaxIter',100000, ...
  'Display','iter','MaxFunEvals',5000); 
%solve using simplex
[param_estimates,fval,exitflag] = fminsearch('simplex_funct',param_start,options)


%function evaluation at starting parameters;
feval_new = simplex_funct(param_start);

%algorithm parameters
alpha = 1; gamma = 1; beta = 0.5; tau = 0.5; s = 1/2;

%form initial simplex
simplex_start = [param_start [param_start(1) + s; param_start(2)] [param_start(1); param_start(2) + s]]

%sort intial simplex by feval values;
for j = 1:(J+1);
feval(j,1) = simplex_funct(simplex_start(:,j));
end;
[feval_sorted,order] = sort(feval);
for j = 1:(J+1)
simplex_new(:,j) = simplex_start(:,order(j));
end;

iter_max = 100;
crit1 = 100;
crit2 = 100;
iter = 0;
param0 = 10*ones(J,1);
while crit1 > 0.001 |  crit2 > 0.001;   
iter = iter + 1;    
if iter > iter_max;
    break;
else
      
%best point
param0 = simplex_new(:,1);
feval0 = simplex_funct(param0);

%next to worst point
paramJ1 = simplex_new(:,J);
fevalJ1 = simplex_funct(paramJ1);

%worst point
paramJ = simplex_new(:,J+1);
fevalJ = simplex_funct(paramJ);

%centroid
M = 1/J * [param0 + paramJ1];

%reflection point
paramR = M + alpha * (M - paramJ);
fevalR = simplex_funct(paramR);

%case 1
if fevalR < feval0; 
display('case 1');;    

%expansion point
paramE = paramR + gamma * (paramR - M);
fevalE = simplex_funct(paramE);

if fevalE < feval0;
simplex_new(:,J+1) = paramE;

else %fevalE >= feval0;
simplex_new(:,J+1) = paramR; 

end;

%case 2
elseif fevalR < fevalJ1;
display('case 2');;       
    
simplex_new(:,J+1) = paramR;


%case 3
else 
 
display('case 3');;  

%contraction point
if fevalR < fevalJ;
    param_tilde = paramR;
else
    param_tilde = paramJ;
end

paramC = M + beta * (param_tilde - M);
fevalC = simplex_funct(paramC);

if fevalC < fevalJ;
    
simplex_new(:,J+1) = paramC;

%shrink
else
simplex_new = simplex_new * tau;
end;   

end; %cases

%sort new simplex by feval values;
simplex_old = simplex_new;
for j = 1:(J+1);
feval(j,1) = simplex_funct(simplex_old(:,j));
end;
[feval_sorted,order] = sort(feval);
for j = 1:(J+1)
simplex_new(:,j) = simplex_old(:,order(j));
end;

crit1 = max(max(abs(simplex_new(:,1) - simplex_new(:,2:J))));
crit2 = max(abs(feval_sorted(1) - feval_sorted(2:J)));

end; %iteration max
end; %crit

%estimates
param_estimates = simplex_new(:,1)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Example 5: Simulated Annealing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;

display('Simulated Annealing');;

%paramaters;
lambda = 1;
T = 0.01;
epsilon = 0.00001;
param_start = 1
param_start = log(param_start); %transform for problem

%maximum number of draws for each iteration
R = 500;

%maximum number of iterations;
iter_max = 1000;

rand('seed',1035);
draws = rand(R,iter_max,2);

z_draws = norminv(draws(:,:,1));
u_draws = draws(:,:,2);

%initial conditions
iter = 0;
r = 1;
crit = epsilon+1;
param_new = param_start;
f_param_old = sim_annealing_funct(param_new);
f_param_new = f_param_old + 10;

%count total objective function evaluations;
count = 0;

while iter < iter_max && crit > epsilon;
    iter = iter+1;
    
    %re-define current parameter as old
    param_old = param_new;
    f_param_old = f_param_new;
    
for r=1:R;
    %create a new parameter
    param_new = param_old + z_draws(r,iter) * lambda;
    f_param_new = sim_annealing_funct(param_new);
    M = (f_param_new - f_param_old) / abs(f_param_old);
    count = count + 1;
    if r == R
        display('max R reached');;
    end;
if f_param_new < f_param_old || M < T * u_draws(r,iter);
break; %accept current parameter of param_new and f_param_new
end; 
end; %r loop

%evaluate stopping criteria
crit = abs(f_param_old - f_param_new) / abs(f_param_old);

end; %while loop

%estimates
param_estimate = exp(param_new)
f_param_new
iter
count



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Example 6: Penalty Function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;

global penalty;

display('Penalty Function');;

%optimization options   
options = optimset('LargeScale','off','TolFun', 0.0001, 'TolX', 0.0001, 'MaxIter',100000, ...
  'Display','iter','MaxFunEvals',1000); 

%penalty function;
penalty = 100;

param_start = 1

display('Without Penalty Function');;

display('case 1');;
%solve using simplex (without penalty function)
[param_estimates,fval,exitflag] = fminsearch('newton_func',param_start,options)

display('case 2');;
%solve using newton-raphson (without penalty function)
[param_estimates,fval,exitflag] = fminunc('newton_func',param_start,options)

display('With Penalty Function');;

display('case 3');;
%solve using simplex (with penalty function)
[param_estimates,fval,exitflag] = fminsearch('newton_func_penalty',param_start,options)

display('case 4');;
%solve using newton-raphson (with penalty function)
[param_estimates,fval,exitflag] = fminunc('newton_func_penalty',param_start,options)


param_start = -1

display('Without Penalty Function');;

display('case 5');;
%solve using simplex (without penalty function)
[param_estimates,fval,exitflag] = fminsearch('newton_func',param_start,options)

display('case 6');;
%solve using newton-raphson (without penalty function)
[param_estimates,fval,exitflag] = fminunc('newton_func',param_start,options)

display('With Penalty Function');;

display('case 7');;
%solve using simplex (with penalty function)
[param_estimates,fval,exitflag] = fminsearch('newton_func_penalty',param_start,options)

display('case 8');;
%solve using newton-raphson (with penalty function)
[param_estimates,fval,exitflag] = fminunc('newton_func_penalty',param_start,options)



diary off;