Update of /cvsroot/octave/octave-forge/main/econometrics
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv15337/main/econometrics
Added Files:
average_moments.m gmm_estimate.m gmm_example.m gmm_obj.m
gmm_results.m gmm_variance_inefficient.m gmm_variance.m
Log Message:
initial commit of gmm functions
--- NEW FILE: gmm_results.m ---
# Copyright (C) 2003,2004 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
function [theta, V, obj_value] = gmm_results(theta, data, weight, moments, momentargs, names, title, unscale, control)
if nargin < 9
[theta, obj_value, convergence] = gmm_estimate(theta, data, weight, moments, momentargs);
else
[theta, obj_value, convergence] = gmm_estimate(theta, data, weight, moments, momentargs, control);
endif
m = feval(moments, theta, data, momentargs); # find out how many obsns. we have
n = rows(m);
if convergence == 1
convergence="Normal convergence";
else
convergence="No convergence";
endif
V = gmm_variance(theta, data, weight, moments, momentargs);
# unscale results if argument has been passed
# this puts coefficients into scale corresponding to the original data
if nargin > 7
if iscell(unscale)
[theta, V] = unscale_parameters(theta, V, unscale);
endif
endif
[theta, V] = delta_method("Parameterize", theta, {data, moments, momentargs}, V);
n = rows(data);
k = rows(theta);
se = sqrt(diag(V));
printf("\n\n******************************************************\n");
disp(title);
printf("\nGMM Estimation Results\n");
printf("BFGS convergence: %s\n", convergence);
printf("\nObjective function value: %f\n", obj_value);
printf("Observations: %d\n", n);
junk = "X^2 test";
df = rows(weight) - rows(theta);
if df > 0
clabels = str2mat("Value","df","p-value");
a = [n*obj_value, df, 1 - chisquare_cdf(n*obj_value, df)];
printf("\n");
prettyprint(a, junk, clabels);
else
disp("\nExactly identified, no spec. test");
end;
# results for parameters
a =[theta, se, theta./se, 2 - 2*normal_cdf(abs(theta ./ se))];
clabels = str2mat("estimate", "st. err", "t-stat", "p-value");
printf("\n");
prettyprint(a, names, clabels);
printf("******************************************************\n");
endfunction
--- NEW FILE: gmm_example.m ---
# Copyright (C) 2003,2004 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# GMM example file, shows initial consistent estimator,
# estimation of efficient weight, and second round
# efficient estimator
# This also shows how do use data scaling - WHICH YOU SHOULD DO!
1;
# the form a user-written moment function should take
function m = mymoments(theta, data, momentargs)
k = momentargs{1}; # use this so that data can hold dep, indeps, and instr
y = data(:,1);
x = data(:,2:k+1);
w = data(:, k+2:columns(data));
lambda = exp(x*theta);
e = y ./ lambda - 1;
m = dmult(e, w);
endfunction
n = 1000;
k = 5;
x = [ones(n,1) rand(n,k-1)];
w = [x, rand(n,1)];
theta = ones(k,1);
lambda = exp(x*theta);
y = randp(lambda);
[xs, scalecoef] = scale_data(x);
# The arguments for gmm_estimate
theta = zeros(k,1);
data = [y xs w];
weight = eye(columns(w));
moments = "mymoments";
momentargs = {k}; # needed to know where x ends and w starts
# additional args for gmm_results
names = str2mat("theta1", "theta2", "theta3", "theta4", "theta5");
title = "Poisson GMM trial";
control = {100,0,1,1};
# initial consistent estimate: only used to get efficient weight matrix, no screen output
[theta, obj_value, convergence] = gmm_estimate(theta, data, weight, moments, momentargs);
# efficient weight matrix
# this method is valid when moments are not autocorrelated
# the user is reponsible to properly estimate the efficient weight
m = feval(moments, theta, data, momentargs);
weight = inverse(cov(m));
# second round efficient estimator
gmm_results(theta, data, weight, moments, momentargs, names, title, scalecoef, control);
--- NEW FILE: gmm_variance.m ---
# Copyright (C) 2003,2004 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# GMM variance, which assumes weights are optimal
function V = gmm_variance(theta, data, weight, moments, momentargs)
D = numgradient("average_moments", {theta, data, moments, momentargs});
D = D';
m = feval(moments, theta, data, momentargs); # find out how many obsns. we have
n = rows(m);
V = (1/n)*inv(D*weight*D');
endfunction
--- NEW FILE: gmm_obj.m ---
# Copyright (C) 2003,2004 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# The GMM objective function
# This is scaled so that it converges to a
# finite number. To get the chi-square specification
# test you need to multiply by n (the sample size)
function obj_value = gmm_obj(theta, data, weight, moments, momentargs)
m = average_moments(theta, data, moments, momentargs);
obj_value = m' * weight *m;
if (((abs(obj_value) == Inf)) || (isnan(obj_value)))
obj_value = realmax;
endif
endfunction
--- NEW FILE: gmm_variance_inefficient.m ---
# Copyright (C) 2003,2004 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# GMM variance, which assumes weights are not optimal
function V = gmm_variance_inefficient(theta, data, weight, omega, moments, momentargs)
D = NumGradient("average_moments", {theta, data, moments, momentargs});
D = D';
m = feval(moments, theta, data, momentargs); # find out how many obsns. we have
n = rows(m);
J = D*weight*D';
J = inv(J);
I = D*weight*omega*weight*D';
V = (1/n)*J*I*J;
endfunction
--- NEW FILE: gmm_estimate.m ---
# Copyright (C) 2003,2004 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# performs the minimization
function [theta, obj_value, convergence] = gmm_estimate(theta, data, weight, moments, momentargs, control)
gmmargs = {theta, data, weight, moments, momentargs};
if nargin < 6
[theta, obj_value, convergence] = bfgsmin("gmm_obj", gmmargs);
else
[theta, obj_value, convergence] = bfgsmin("gmm_obj", gmmargs, control);
endif
endfunction
--- NEW FILE: average_moments.m ---
# Copyright (C) 2003,2004 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# average moments (separate function so it can be differentiated)
function m = average_moments(theta, data, moments, momentargs)
m = feval(moments, theta, data, momentargs);
m = mean(m)'; # returns Gx1 moment vector
endfunction

Update of /cvsroot/octave/octave-forge/main/econometrics
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv30326/main/econometrics
Added Files:
delta_method.m mle_estimate.m mle_example.m mle_obj.m
mle_results.m mle_variance.m parameterize.m prettyprint_c.m
prettyprint.m scale_data.m unscale_parameters.m
Log Message:
initial commit of econometrics functions
--- NEW FILE: scale_data.m ---
# Copyright (C) 2003,2004 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Standardizes and normalizes data matrix,
# primarily for use by BFGS
function [zz, scalecoefs] = scale_data(z);
n = rows(z);
k = columns(z);
# Scale data
s = std(z)';
test = s != 0;
s = s + (1 - test); # don't scale if column is a constant (avoid div by zero)
A = diag(1 ./ s);
# De-mean all variables except constant, if a constant is present
test = std(z(:,1)) != 0;
if test
warning("ScaleData: no constant present - only scale - no de-mean");
bb = zeros(n,k);
b = zeros(k,1);
else
b = -mean(z)';
test = std(z)' != 0;
# don't take out mean if the column is a constant, to preserve identification
b = b .* test;
b = A*b;
bb = dmult(b, ones(k,n))';
endif
zz = z*A + bb;
scalecoefs = {A,b};
endfunction
--- NEW FILE: prettyprint.m ---
# Copyright (C) 2003,2004 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# this prints matrices with row and column labels
function prettyprint(mat, rlabels, clabels)
# left pad the column labels
a = columns(rlabels);
for i = 1:a
printf(" ");
endfor
printf(" ");
# print the column labels
clabels = [" ";clabels]; # pad to 8 characters wide
clabels = strjust(clabels,"right");
k = columns(mat);
for i = 1:k
printf("%s ",clabels(i+1,:));
endfor
# now print the row labels and rows
printf("\n");
k = rows(mat);
for i = 1:k
if isstr(rlabels(i,:))
printf(rlabels(i,:));
else
printf("%i", rlabels(i,:));
endif
printf(" %10.3f", mat(i,:));
printf("\n");
endfor
endfunction
--- NEW FILE: mle_example.m ---
# Copyright (C) 2003,2004 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Example to show how to use MLE functions
# Example likelihood function, no score
function [log_density, score] = poisson_no_score(theta, data, otherargs)
y = data(:,1);
x = data(:,2:columns(data));
lambda = exp(x*theta);
log_density = -lambda + y .* (x*theta) - lgamma(y+1);
score = "na";
endfunction
# Example likelihood function, with score
function [log_density, score] = poisson_with_score(theta, data, otherargs)
y = data(:,1);
x = data(:,2:columns(data));
lambda = exp(x*theta);
log_density = -lambda + y .* (x*theta) - lgamma(y+1);
score = dmult(y - lambda,x);
endfunction
# Generate data
n = 1000; # how many observations?
# the explanatory variables: note that they have unequal scales
x = [ones(n,1) rand(n,1) randn(n,1)];
theta = 1:3; # true coefficients are 1,2,3
theta = theta';
lambda = exp(x*theta);
y = randp(lambda); # generate the dependent variable
####################################
# define arguments for mle_results #
####################################
# starting values
theta = zeros(3,1);
# data
data = [y, x];
# name of model to estimate
model = "poisson_with_score";
# placeholder, poisson model has no additional args
modelargs = {};
# parameter names
names = str2mat("beta1", "beta2", "beta3");
title = "Poisson MLE trial"; # title for the run
# controls for bfgsmin: 30 iterations is not always enough for convergence
control = {30,0,1,1};
# This displays the results
printf("\n\nanalytic score, unscaled data\n");
[theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, title, 0, control);
# This just calculates the results, no screen display
printf("\n\nanalytic score, unscaled data, no screen display\n");
theta = zeros(3,1);
[theta, obj_value, convergence] = mle_estimate(theta, data, model, modelargs, control);
printf("obj_value = %f, to confirm it worked\n", obj_value);
# This show how to scale data during estimation, but get results
# for data in original units (recommended to avoid conditioning problems)
# This usually converges faster, depending upon the data
printf("\n\nanalytic score, scaled data\n");
[scaled_x, unscale] = scale_data(x);
data = [y, scaled_x];
theta = zeros(3,1);
[theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, title, unscale, control);
# Example using numeric score
printf("\n\nnumeric score, scaled data\n");
theta = zeros(3,1);
model = "poisson_no_score";
[theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, title, unscale, control);
--- NEW FILE: mle_results.m ---
# Copyright (C) 2003,2004 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# report results
function [theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, title, unscale, control)
if nargin < 8
[theta, obj_value, convergence] = mle_estimate(theta, data, model, modelargs);
else
[theta, obj_value, convergence] = mle_estimate(theta, data, model, modelargs, control);
endif
V = mle_variance(theta, data, model, modelargs);
# unscale results if argument has been passed
# this puts coefficients into scale corresponding to the original modelargs
if (nargin > 5)
if iscell(unscale) # don't try it if unscale is simply a placeholder
[theta, V] = unscale_parameters(theta, V, unscale);
endif
endif
[theta, V] = delta_method("parameterize", theta, {data, model, modelargs}, V);
n = rows(data);
k = rows(V);
se = sqrt(diag(V));
if convergence == 1
convergence="Normal convergence";
elseif convergence == 2
convergence="No convergence";
elseif convergence == -1
convergence = "Max. iters. exceeded";
endif
printf("\n\n******************************************************\n");
disp(title);
printf("\nMLE Estimation Results\n");
printf("BFGS convergence: %s\n\n", convergence);
printf("Average Log-L: %f\n", obj_value);
printf("Observations: %d\n", n);
a =[theta, se, theta./se, 2 - 2*normal_cdf(abs(theta ./ se))];
clabels = str2mat("estimate", "st. err", "t-stat", "p-value");
printf("\n");
if names !=0
prettyprint(a, names, clabels);
else prettyprint_c(a, clabels);
endif
printf("\nInformation Criteria \n");
caic = -2*n*obj_value + rows(theta)*(log(n)+1);
bic = -2*n*obj_value + rows(theta)*log(n);
aic = -2*n*obj_value + 2*rows(theta);
infocrit = [caic, bic, aic];
printf("CAIC : %8.4f Avg. CAIC: %8.4f\n", caic, caic/n);
printf(" BIC : %8.4f Avg. BIC: %8.4f\n", bic, bic/n);
printf(" AIC : %8.4f Avg. AIC: %8.4f\n", aic, aic/n);
printf("******************************************************\n");
endfunction
--- NEW FILE: delta_method.m ---
# Copyright (C) 2003,2004 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Computes Delta method mean and covariance of a nonlinear
# transformation defined by "func"
function [theta_transf, vartheta] = delta_method(func, theta, otherargs, vartheta)
theta_transf = feval(func, theta, otherargs);
D = numgradient(func, {theta, otherargs});
var_transf = D * vartheta * D';
endfunction
--- NEW FILE: unscale_parameters.m ---
# Copyright (C) 2003,2004 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Unscales parameters that were estimated using scaled data
# primarily for use by BFGS
function [theta_us, vartheta_us] = unscale_parameters(theta, vartheta, scalecoefs);
k = rows(theta);
A = nth(scalecoefs, 1);
b = nth(scalecoefs, 2);
kk = rows(b);
B = zeros(kk-1,kk);
B = [b'; B];
C = A + B;
# allow for parameters that aren't associated with x
if (k > kk)
D = zeros(kk, k - kk);
C = [C, D; D', eye(k - kk)];
endif;
theta_us = C*theta;
vartheta_us = C * vartheta * C';
endfunction
--- NEW FILE: mle_obj.m ---
# Copyright (C) 2003,2004 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# this takes a general model and calculates average log likelihood
function [obj_value, score] = mle_obj(theta, data, model, modelargs)
[obj_value, score] = feval(model, theta, data, modelargs);
obj_value = - mean(obj_value);
# let's bullet-proof this in case the model goes nuts
if (((abs(obj_value) == Inf)) || (isnan(obj_value)))
obj_value = realmax;
endif
if isnumeric(score) score = - mean(score)'; endif
endfunction
--- NEW FILE: parameterize.m ---
# Copyright (C) 2003,2004 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
#
# This is an empty function, provided so that
# delta_method will work as is. Replace it with
# the parameter transformations your models use.
# Note: you can let "otherargs" contain the model
# name so that this function can do parameterizations
# for a variety of models
function theta = parameterize(theta, otherargs)
endfunction
--- NEW FILE: mle_estimate.m ---
# Copyright (C) 2003,2004 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# call the minimizing routine
function [theta, obj_value, convergence] = mle_estimate(theta, data, model, modelargs, control)
if nargin == 3
[theta, obj_value, convergence] = bfgsmin("mle_obj", {theta, data, model, modelargs});
else
[theta, obj_value, convergence] = bfgsmin("mle_obj", {theta, data, model, modelargs}, control);
endif
obj_value = - obj_value; # recover from minimization rather than maximization
endfunction
--- NEW FILE: prettyprint_c.m ---
# Copyright (C) 2003,2004 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# this prints matrices with column labels but no row labels
function prettyprint_c(mat, clabels)
printf(" ");
# print the column labels
clabels = [" ";clabels]; # pad to 8 characters wide
clabels = strjust(clabels,"right");
k = columns(mat);
for i = 1:k
printf("%s ",clabels(i+1,:));
endfor
# now print the row labels and rows
printf("\n");
k = rows(mat);
for i = 1:k
printf(" %8.3f", mat(i,:));
printf("\n");
endfor
endfunction
--- NEW FILE: mle_variance.m ---
# Copyright (C) 2003,2004 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# Copyright (C) 2003 Michael Creel michael.creel@...
# under the terms of the GNU General Public License.
# The GPL license is in the file COPYING
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
# sandwich form of var-cov matrix
function [V,scorecontribs,J_inv] = mle_variance(theta, data, model, modelargs)
scorecontribs = numgradient(model, {theta, data, modelargs});
n = rows(scorecontribs);
I = scorecontribs'*scorecontribs / n;
J = numhessian("mle_obj", {theta, data, model, modelargs});
J_inv = inverse(J);
V = J_inv*I*J_inv/n;
endfunction

Community

Help

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

CountryState

JavaScript is required for this form.

I agree to receive quotes, newsletters and other information from sourceforge.net and its partners regarding IT services and products. I understand that I can withdraw my consent at any time. Please refer to our Privacy Policy or Contact Us for more details