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

clc
clear

global x T MA

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

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

T = 50;                                                                     % Number of periods.
MA = 2;                                                                     % MA order.
sigma = 2;                                                                  % Standard deviation.
beta1 = 0.5;                                                                % MA coefficient.

if MA == 1                                                                  % For MA(1) process.
    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.
    
elseif MA == 2                                                              % For MA(2) process.
    beta2 = 0.4;                                                            % MA coefficient.
    sim = 5;                                                                % Number of simulations.
    results_mat = zeros(sim,4);                                             % Results matrix.
    se_mat = zeros(sim,4);                                                  % Standard errors.
    t_mat = zeros(sim,4);                                                   % t-stats.
    
end

% Begin Simulation.

for s = 1:sim
    
    % Generate the data using the model.
    %               MA(1):   x(t) = u(t) + beta1*u(t-1)
    %               MA(2):   x(t) = u(t) + beta1*u(t-1) + beta2*u(t-2)   
    
    u = normrnd(0,sigma,T,1);
    x = zeros(T,1);
    
    if MA == 1                                                              % For MA(1) process.
        
        x(1) = u(1) + beta1*normrnd(0,sigma,1,1);
        
        for j = 2:T
            
            x(j) = u(j) + beta1*u(j-1);
            
        end
        
        init = [1 1 0.1];                                                   % Initial values.
        
    elseif MA == 2                                                          % For MA(2) process.
        
        c = normrnd(0,sigma,1,1);
        x(1) = u(1) + beta1*c + beta2*normrnd(0,sigma,1,1);
        x(2) = u(2) + beta1*u(1) + beta2*c;
        
        for j = 3:T
            
            x(j) = u(j) + beta1*u(j-1) + beta2*u(j-2);
            
        end
        
        init = [1 1 0.1 0.1];                                               % Initial values.
        
    end
    
    options = optimset('Display','off');                                    % Turn off the display.
    [b_mle,~,~,~,~,hess] = fminunc('logl_MA',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))
    if MA == 2
        fprintf('beta2   %0.4f       %0.4f         %0.4f    \n', results_mat(k,4), se_mat(k,4), t_mat(k,4))           % Extra line for MA(2).
    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 MA == 2
fprintf('beta2   %0.4f       %0.4f     \n', mean(results_mat(:,4)), std(results_mat(:,4)))          % Extra line for MA(2).
end
fprintf('\n')
