%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Econometrics 1
% Problem Set 11
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This code simulates and estimates the parameters of an AR(1).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc
clear

global x T AR

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

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

sim = 2;                                                                   % Number of simulations.

if AR == 1
    
    results_mat = zeros(sim,3);                                             % Results matrix.
    se_mat = zeros(sim,3);                                                  % Standard errors.
    t_mat = zeros(sim,3);                                                   % t-stats.
    se_opg = zeros(sim,3);                                                  % OPG Standard errors.
    
elseif AR == 2
    
    results_mat = zeros(sim,4);                                             % Results matrix.
    se_mat = zeros(sim,4);                                                  % Standard errors.
    t_mat = zeros(sim,4);                                                   % t-stats.
    beta2 = 0.4;                                                            % AR coefficient.
    se_opg = zeros(sim,4);                                                  % OPG Standard errors.
    
end

% Begin Simulation.

for s = 1:sim
    
    % Generate the data using the model.
    %              AR(1):  x(t) = mu + beta1*x(t-1) + u(t)
    %              AR(2):  x(t) = mu + beta1*x(t-1) + beta(2)*x(t-2) + u(t)
    
    u = normrnd(0,sigma,T,1);
    x = zeros(T,1);
    
    if AR == 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 = [0 1 0];                                                     % Initial values.
        
    elseif AR == 2
    
        x(1) = mu/(1-beta1-beta2) + u(1)*sqrt((1-beta2)...
                        /(1-beta2-beta1^2*(1-beta2)+beta2^2*(1-beta2)));  
        x(2) = mu/(1-beta1-beta2) + beta1/(1-beta2)*(x(1)...
                    -(mu/(1-beta1-beta2)))+sqrt((1-beta2)...
                       /(1-beta2-beta1^2*(1-beta2)+beta2^2*(1-beta2)))*...
                            sqrt(1-(beta1/(1-beta2))^2)*u(2);
       
       %x(1) = mu + u(1);                                                   % Head-ache free first terms.
       %x(2) = mu + beta1*x(1) + u(2);
       
            for j = 3:T
            
                x(j) = mu + beta1*x(j-1) + beta2*x(j-2) + u(j);                                
            
            end
        
        init = [1 0.5 0.4 0.1];                                             % Initial values.
    
    end
           
    options = optimset('Display','off','algorithm','quasi-newton');         % Turn off the display.
    [b_mle,~,~,~,~,hess] = fminunc('logl_AR',init,options);                 % Minimization.
    
    % OPG estimator. This takes a while.
    
    if AR == 1
        B = sym('b0',[1 2 3]);                                              % Declare symbolic variables.
    elseif AR == 2
        B = sym('b0',[1 2 3 4]);                                            
    end
    
    sL = logl_AR(B);                                                        % Evaluate the loglikelihood at the symbolic variables.
    sL = symfun(sL,symvar(sL));                                             % Convert the Matlab function into a symbolic function.
    sL = children(sL);                                                      % Break up the loglikelihood function into its contributions.
    G = jacobian(sL);                                                       % Matrix of contributions to the loglikelihood function
    G = matlabFunction(G);                                                  % Convert the matrix of contributions back into a Matlab function.
    G = G(b_mle(1),b_mle(2),b_mle(3),b_mle(4));                             % Evaluate at the MLE estimates.
    I = G'*G;                                                               % Information Matrix.
    OPG = inv(I);                                                           % OPG estimator of covariance matrix.
    
    results_mat(s,:) = b_mle';                                              % Store estimates.
    se_mat(s,:) = sqrt(abs(diag(inv(hess))))';                              % Store standard errors.
    t_mat(s,:) = results_mat(s,:)./se_mat(s,:);                             % Store t-stats.
    se_opg(s,:) = sqrt(abs(diag(OPG)))';                                    % Store OPG standard errors.
    
end

% Display the results of each simulation.

for k = 1: size(results_mat,1)
    fprintf('Simulation %d \n',k)
    fprintf('          b            SE             t            OPG \n')
    fprintf('mu      %0.4f       %0.4f         %0.4f         %0.4f\n', results_mat(k,1), se_mat(k,1), t_mat(k,1), se_opg(k,1))
    fprintf('sigma   %0.4f       %0.4f         %0.4f         %0.4f\n', results_mat(k,2), se_mat(k,2), t_mat(k,2), se_opg(k,2))
    fprintf('beta1   %0.4f       %0.4f         %0.4f         %0.4f\n', results_mat(k,3), se_mat(k,3), t_mat(k,3), se_opg(k,3))
    if AR == 2
        fprintf('beta2   %0.4f       %0.4f         %0.4f         %0.4f\n', results_mat(k,4), se_mat(k,4), t_mat(k,4), se_opg(k,4))
    end
    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)))
    if AR == 2
        fprintf('beta2   %0.4f       %0.4f         %0.4f    \n', mean(results_mat(:,4)), std(results_mat(:,4)))
    end
fprintf('\n')

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

%Add OLS here
