%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Xavier Martin G. Bautista
% Econometrics 2
% Problem Set 3
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This code simulates and estimates the parameters of an AR(1).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc
clear

global x T

fprintf('--------------------------------------------------------------\n')
fprintf('Name: Xavier Martin G. Bautista \n')
fprintf('Econometrics 2 \n')
fprintf('Problem Set 3 \n')
fprintf('--------------------------------------------------------------\n')
fprintf('\n')

% Set up the true parameters and placeholders for the results.

T = 50;                                                                     % Number of periods.
beta1 = 0.5;                                                                % AR coefficient.
sigma = 2;                                                                  % Constant.
mu = 3;

sim = 10;                                                                    % Number of simulations.
results_mat = zeros(sim,3);                                                 % Results matrix.
se_mat = zeros(sim,3);                                                      % Standard errors.
t_mat = zeros(sim,3);                                                       % t-stats.

% Begin Simulation.

for s = 1:sim
    
    % Generate the data using the model.
    %               AR(1):   x(t) = mu + beta1*x(t-1) + u(t)
    
    u = normrnd(0,sigma,T,1);
    x = zeros(T,1);
    
    x(1) = mu/(1-beta1) + u(1)/sqrt(1-beta1^2);                               
    
        for j = 2:T
            
            x(j) = mu + beta1*x(j-1) + u(j);                                
            
        end
        
    init = [1 1 0.1];                                                       % Initial values.
           
    options = optimset('Display','off');                                    % Turn off the display.
    [b_mle,~,~,~,grad,hess] = fminunc('logl_AR',init,options);              % Minimization.
    
    results_mat(s,:) = b_mle';                                              % Store estimates.
    se_mat(s,:) = sqrt(diag(inv(hess)))';                                   % Store standard errors.
    t_mat(s,:) = results_mat(s,:)./se_mat(s,:);                             % Store t-stats.
    
end

% Display the results of each simulation.

for k = 1: size(results_mat,1)
    fprintf('Simulation %d \n',k)
    fprintf('          b            SE             t \n')
    fprintf('mu      %0.4f       %0.4f         %0.4f    \n', results_mat(k,1), se_mat(k,1), t_mat(k,1))
    fprintf('sigma   %0.4f       %0.4f         %0.4f    \n', results_mat(k,2), se_mat(k,2), t_mat(k,2))
    fprintf('beta1   %0.4f       %0.4f         %0.4f    \n', results_mat(k,3), se_mat(k,3), t_mat(k,3))
    fprintf('\n')
end

fprintf('--------------------------------------------------------------\n')
fprintf('\n')

% Display the empirical means and standard deviations.

fprintf('Empirical Results \n')
fprintf('         Mean       Std.Dev   \n')
fprintf('mu      %0.4f       %0.4f     \n', mean(results_mat(:,1)), std(results_mat(:,1)))
fprintf('sigma   %0.4f       %0.4f     \n', mean(results_mat(:,2)), std(results_mat(:,2)))
fprintf('beta1   %0.4f       %0.4f     \n', mean(results_mat(:,3)), std(results_mat(:,3)))
fprintf('\n')