Difference between revisions of "MaxLikCode"

From ECLR
Jump to: navigation, search
(Created page with "== Maximum Likelihood codes == = MLse = This code simulates a linear model, estimates it by PLS and ML and calculates standard errors <source> % Code to % a) simulate line...")
 
(Maximum Likelihood codes)
Line 1: Line 1:
== Maximum Likelihood codes ==  
+
== Maximum Likelihood codes ==
 +
This Section contains a number of codes that are used in the Maximum Likelihood Section.
  
 
= MLse =
 
= MLse =

Revision as of 10:32, 18 October 2012

Maximum Likelihood codes

This Section contains a number of codes that are used in the Maximum Likelihood Section.

MLse

This code simulates a linear model, estimates it by PLS and ML and calculates standard errors

% Code to
% a) simulate linear model, y_i = 0.2 + 0.6 x_1i -1.5 x_2i + err_i
%       x_1i and x_s1 come from N(0,1)
%       err_i ~ N(0,sd^2)
% b) estimate it by OLS
% c) estimate ot by ML
% d) estimate different ML standard errors
%
% This code requires the following Functions:
% OLSest, nll_lin, HessMp, gradp
clc
clear all

T   = 1000;                     % set sample size
b0  = [0.2; 0.6; -1.5];         % set parameter values
sd  = 0.5;                      % set error standard deviation
x   = [ones(T,1) randn(T,2)];   % simulate X matrix
err = randn(T,1)*sd;            % simulate error terms
y = x*b0 + err;                 % calculate y_i s

[b,bse,res,n,rss,r2] = OLSest(y,x,1);   % OLS estimation

datamat = [y x];                        % define data matrix for use in nll_lin
theta0 = [mean(y); 0; 0; std(y)];       % this sets the initial parameter vector
options = optimset; % sets optimisation options to default
[betaopt] = fminunc(@nll_lin,theta0,options,datamat,1);
 
H = HessMp(@nll_lin,betaopt,datamat,1); % this returns the negative of the Hessian
g = gradp(@nll_lin,betaopt,datamat,0);  % this returns a (T x size(theta0,1)) matrix of gradients
J = (g'*g);                             % calculates the OPG

se_H = sqrt(diag(inv(H)));
se_J = sqrt(diag(inv(J)));
se_SW = sqrt(diag(inv(H*inv(J)*H)));    % Sandwich variance covariance

disp('   Est     se(OLS)     se(H)     se(J)     se(SW)');
disp([betaopt [bse;0] se_H se_J se_SW]);