octave_packages/econometrics-1.0.8/mle_example.m

   1 ## Copyright (C) 2003,2004  Michael Creel <michael.creel@uab.es>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 2 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
  15
  16 ## Example to show how to use MLE functions
  17
  18
  19
  20 # Generate data
  21 n = 1000; # how many observations?
  22
  23 # the explanatory variables: note that they have unequal scales
  24 x = [ones(n,1) -rand(n,1) randn(n,1)];
  25 theta = 1:3; # true coefficients are 1,2,3
  26 theta = theta';
  27
  28 lambda = exp(x*theta);
  29 y = poissrnd(lambda); # generate the dependent variable
  30
  31 ####################################
  32 # define arguments for mle_results #
  33 ####################################
  34
  35 # starting values
  36 theta = zeros(3,1);
  37 # data
  38 data = [y, x];
  39 # name of model to estimate
  40 model = "poisson";
  41 modelargs = {0}; # if this is zero the function gives analytic score, otherwise not
  42 # parameter names
  43 names = char("beta1", "beta2", "beta3");
  44 mletitle = "Poisson MLE trial"; # title for the run
  45
  46 # controls for bfgsmin: 30 iterations is not always enough for convergence
  47 control = {50,0};
  48
  49 # This displays the results
  50 printf("\n\nanalytic score, unscaled data\n");
  51 [theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, mletitle, 0, control);
  52
  53 # This just calculates the results, no screen display
  54 printf("\n\nanalytic score, unscaled data, no screen display\n");
  55 theta = zeros(3,1);
  56 [theta, obj_value, convergence] = mle_estimate(theta, data, model, modelargs, control);
  57 printf("obj_value = %f, to confirm it worked\n", obj_value);
  58
  59 # This show how to scale data during estimation, but get results
  60 # for data in original units (recommended to avoid conditioning problems)
  61 # This usually converges faster, depending upon the data
  62 printf("\n\nanalytic score, scaled data\n");
  63 [scaled_x, unscale] = scale_data(x);
  64 data = [y, scaled_x];
  65 theta = zeros(3,1);
  66 [theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, mletitle, unscale, control);
  67
  68 # Example using numeric score
  69 printf("\n\nnumeric score, scaled data\n");
  70 theta = zeros(3,1);
  71 modelargs = {1}; # set the switch for no score
  72 [theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, mletitle, unscale, control);
  73
  74 # Example doing estimation in parallel on a cluster (requires MPITB)
  75 # uncomment the following if you have MPITB installed
  76 # theta = zeros(3,1);
  77 # nslaves = 1;
  78 # title = "MLE estimation done in parallel";
  79 # [theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, mletitle, unscale, control, nslaves);