Menu
▾
▴
[Octave-cvsupdate] octave-forge/main/econometrics average_moments.m,1.4,1.5 gmm_estimate.m,1.4,1.5 gmm_example.m,1.5,1.6 gmm_results.m,1.5,1.6 mle_example.m,1.5,1.6 mle_obj.m,1.4,1.5 mle_obj_nodes.m,1.2,1.3 mle_results.m,1.4,1.5 mle_variance.m,1.3,1.4 nls_obj.m,1.2,1.3 nls_obj_nodes.m,1.2,1.3 sum_moments_nodes.m,1.2,1.3
|
From: Michael C. <mc...@us...> - 2005-05-30 15:16:23
|
Update of /cvsroot/octave/octave-forge/main/econometrics In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv10062/main/econometrics Modified Files: average_moments.m gmm_estimate.m gmm_example.m gmm_results.m mle_example.m mle_obj.m mle_obj_nodes.m mle_results.m mle_variance.m nls_obj.m nls_obj_nodes.m sum_moments_nodes.m Log Message: cleanups and documentation Index: gmm_results.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/econometrics/gmm_results.m,v retrieving revision 1.5 retrieving revision 1.6 diff -u -d -r1.5 -r1.6 --- gmm_results.m 25 May 2005 13:04:33 -0000 1.5 +++ gmm_results.m 30 May 2005 15:16:11 -0000 1.6 @@ -18,24 +18,30 @@ # gmm_results(theta, data, weight, moments, momentargs, names, title, unscale, control, nslaves) # # inputs: -# theta: column vector initial parameters -# data: data matrix -# weight: the GMM weight matrix -# moments: name of function computes the moments (should return nXg matrix of contributions) -# momentargs: (cell) additional inputs needed to compute moments. May be empty ("") -# names: vector of parameter names, e.g., use names = str2mat("param1", "param2"); -# title: string, describes model estimated -# unscale: (optional) cell that holds means and std. dev. of data (see scale_data) -# control: (optional) BFGS or SA controls (see bfgsmin and samin). May be empty (""). -# nslaves: (optional) number of slaves if executed in parallel (requires MPITB) +# theta: column vector initial parameters +# data: data matrix +# weight: the GMM weight matrix +# moments: name of function computes the moments +# (should return nXg matrix of contributions) +# momentargs: (cell) additional inputs needed to compute moments. +# May be empty ("") +# names: vector of parameter names +# e.g., names = str2mat("param1", "param2"); +# title: string, describes model estimated +# unscale: (optional) cell that holds means and std. dev. of data +# (see scale_data) +# control: (optional) BFGS or SA controls (see bfgsmin and samin). May be empty (""). +# nslaves: (optional) number of slaves if executed in parallel +# (requires MPITB) # # outputs: # theta: GMM estimated parameters -# V: estimate of covariance of parameters. Assumes the weight matrix is optimal +# V: estimate of covariance of parameters. Assumes the weight matrix +# is optimal (inverse of covariance of moments) # obj_value: the value of the GMM objective function # -# Please see gmm_example for information on how to use this - +# please type "gmm_example" while in octave to see an example + function [theta, V, obj_value] = gmm_results(theta, data, weight, moments, momentargs, names, title, unscale, control, nslaves) Index: nls_obj.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/econometrics/nls_obj.m,v retrieving revision 1.2 retrieving revision 1.3 diff -u -d -r1.2 -r1.3 --- nls_obj.m 25 May 2005 13:04:33 -0000 1.2 +++ nls_obj.m 30 May 2005 15:16:12 -0000 1.3 @@ -24,40 +24,38 @@ n = rows(data); - if nslaves > 0 - global NEWORLD NSLAVES TAG + if nslaves > 0 + global NEWORLD NSLAVES TAG nn = floor(n/(NSLAVES + 1)); # number of obsns per slave # The command that the slave nodes will execute - cmd=['contrib = nls_obj_nodes(theta, data, model, modelargs, nn); ',... - 'MPI_Send(contrib,0,TAG,NEWORLD);']; + cmd=['contrib = nls_obj_nodes(theta, data, model, modelargs, nn); ',... + 'MPI_Send(contrib,0,TAG,NEWORLD);']; # send items to slaves NumCmds_Send({"theta", "nn", "cmd"}, {theta, nn, cmd}); # evaluate last block on master while slaves are busy - obj_value = nls_obj_nodes(theta, data, model, modelargs, nn); + obj_value = nls_obj_nodes(theta, data, model, modelargs, nn); # collect slaves' results contrib = 0.0; # must be initialized to use MPI_Recv - for i = 1:NSLAVES - MPI_Recv(contrib,i,TAG,NEWORLD); + for i = 1:NSLAVES + MPI_Recv(contrib,i,TAG,NEWORLD); obj_value = obj_value + contrib; endfor # compute the average - obj_value = obj_value / n; - score = "na"; # fix this later to allow analytic score in parallel + obj_value = obj_value / n; + score = "na"; # fix this later to allow analytic score in parallel - else # serial version - [contribs, score] = feval(model, theta, data, modelargs); + else # serial version + [contribs, score] = feval(model, theta, data, modelargs); obj_value = mean(contribs); - if isnumeric(score) score = mean(score)'; endif # model passes "na" when score not available - endif + if isnumeric(score) score = mean(score)'; endif # model passes "na" when score not available + endif # let's bullet-proof this in case the model goes nuts - if (((abs(obj_value) == Inf)) || (isnan(obj_value))) - obj_value = realmax/10; - endif + if (((abs(obj_value) == Inf)) || (isnan(obj_value))) obj_value = realmax; endif endfunction Index: mle_obj_nodes.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/econometrics/mle_obj_nodes.m,v retrieving revision 1.2 retrieving revision 1.3 diff -u -d -r1.2 -r1.3 --- mle_obj_nodes.m 25 May 2005 13:04:33 -0000 1.2 +++ mle_obj_nodes.m 30 May 2005 15:16:12 -0000 1.3 @@ -15,7 +15,7 @@ ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA function contrib = mle_obj_nodes(theta, data, model, modelargs, nn) - global NEWORLD NSLAVES + global NEWORLD NSLAVES # Who am I? [info, rank] = MPI_Comm_rank(NEWORLD); @@ -28,7 +28,7 @@ endif data = data(startblock:endblock,:); - contrib = feval(model, theta, data, modelargs); + contrib = feval(model, theta, data, modelargs); contrib = sum(contrib); endfunction Index: mle_example.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/econometrics/mle_example.m,v retrieving revision 1.5 retrieving revision 1.6 diff -u -d -r1.5 -r1.6 --- mle_example.m 25 May 2005 13:04:33 -0000 1.5 +++ mle_example.m 30 May 2005 15:16:11 -0000 1.6 @@ -73,15 +73,8 @@ [theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, title, unscale, control); # Example doing estimation in parallel on a cluster (requires MPITB) -if exist("MPI_Init") - -printf("to do estimation in parallel, you need to have MPITB installed\n... -and your computer must be lambooted. If this is not the case press... - CRTL-C to abort. Pausing 10 seconds\n"); - pause(10); - theta = zeros(3,1); - nslaves = 1; - title = "MLE estimation done in parallel"; - [theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, title, unscale, control, nslaves); -else printf("sorry, MPITB is not installed, can't do estimation in parallel\n"); -endif +# uncomment the following if you have MPITB installed +# theta = zeros(3,1); +# nslaves = 1; +# title = "MLE estimation done in parallel"; +# [theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, title, unscale, control, nslaves); Index: mle_results.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/econometrics/mle_results.m,v retrieving revision 1.4 retrieving revision 1.5 diff -u -d -r1.4 -r1.5 --- mle_results.m 25 May 2005 13:04:33 -0000 1.4 +++ mle_results.m 30 May 2005 15:16:12 -0000 1.5 @@ -40,58 +40,54 @@ # report results function [theta, V, obj_value, infocrit] = mle_results(theta, data, model, modelargs, names, title, unscale, control, nslaves) - if nargin < 9 nslaves = 0; endif # serial by default + if nargin < 9 nslaves = 0; endif # serial by default if nargin < 8 - [theta, obj_value, convergence] = mle_estimate(theta, data, model, modelargs, "", nslaves); - else - [theta, obj_value, convergence] = mle_estimate(theta, data, model, modelargs, control, nslaves); - endif - V = mle_variance(theta, data, model, modelargs); + [theta, obj_value, convergence] = mle_estimate(theta, data, model, modelargs, "", nslaves); + else + [theta, obj_value, convergence] = mle_estimate(theta, data, model, modelargs, control, nslaves); + 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 + # 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); + [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))]; + 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); - clabels = str2mat("estimate", "st. err", "t-stat", "p-value"); - - printf("\n"); - if names !=0 - prettyprint(a, names, clabels); + 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"); + + 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 Index: average_moments.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/econometrics/average_moments.m,v retrieving revision 1.4 retrieving revision 1.5 diff -u -d -r1.4 -r1.5 --- average_moments.m 25 May 2005 13:04:33 -0000 1.4 +++ average_moments.m 30 May 2005 15:16:11 -0000 1.5 @@ -20,30 +20,30 @@ # average moments (separate function so it can be differentiated) function m = average_moments(theta, data, moments, momentargs) - global NSLAVES PARALLEL NEWORLD NSLAVES TAG; + global NSLAVES PARALLEL NEWORLD NSLAVES TAG; n = rows(data); - if PARALLEL - nn = floor(n/(NSLAVES + 1)); # save some work for master + if PARALLEL + nn = floor(n/(NSLAVES + 1)); # save some work for master # The command that the slave nodes will execute - cmd=['contrib = sum_moments_nodes(theta, data, moments, momentargs, nn); ',... - 'MPI_Send(contrib,0,TAG,NEWORLD);']; + cmd=['contrib = sum_moments_nodes(theta, data, moments, momentargs, nn); ',... + 'MPI_Send(contrib,0,TAG,NEWORLD);']; # send items to slaves NumCmds_Send({"theta", "nn", "cmd"},{theta, nn, cmd}); # evaluate last block on master while slaves are busy - m = feval("sum_moments_nodes", theta, data, moments, momentargs, nn); + m = feval("sum_moments_nodes", theta, data, moments, momentargs, nn); # collect slaves' results contrib = zeros(1,columns(m)); - for i = 1:NSLAVES - MPI_Recv(contrib,i,TAG,NEWORLD); + for i = 1:NSLAVES + MPI_Recv(contrib,i,TAG,NEWORLD); m = m + contrib; - endfor - + endfor + m = m'; # we want a column vector, please m = m/n; # average please, not sum @@ -52,5 +52,4 @@ m = mean(m)'; # returns Gx1 moment vector endif - endfunction Index: gmm_example.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/econometrics/gmm_example.m,v retrieving revision 1.5 retrieving revision 1.6 diff -u -d -r1.5 -r1.6 --- gmm_example.m 25 May 2005 13:04:33 -0000 1.5 +++ gmm_example.m 30 May 2005 15:16:11 -0000 1.6 @@ -57,18 +57,9 @@ # Example doing estimation in parallel on a cluster (requires MPITB) -if exist("MPI_Init") - -printf("to do estimation in parallel, you need to have MPITB installed\n... -and your computer must be lambooted. If this is not the case press... - CRTL-C to abort. Pausing 10 seconds\n"); - pause(10); - - nslaves = 1; - theta = zeros(k,1); - nslaves = 1; - title = "GMM estimation done in parallel"; - gmm_results(theta, data, weight, moments, momentargs, names, title, scalecoef, control, nslaves); - -else printf("sorry, MPITB is not installed, can't do estimation in parallel\n"); -endif +# uncomment the following if you have MPITB installed +# nslaves = 1; +# theta = zeros(k,1); +# nslaves = 1; +# title = "GMM estimation done in parallel"; +# gmm_results(theta, data, weight, moments, momentargs, names, title, scalecoef, control, nslaves); Index: mle_obj.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/econometrics/mle_obj.m,v retrieving revision 1.4 retrieving revision 1.5 diff -u -d -r1.4 -r1.5 --- mle_obj.m 25 May 2005 13:04:33 -0000 1.4 +++ mle_obj.m 30 May 2005 15:16:11 -0000 1.5 @@ -26,41 +26,40 @@ if nargin < 5 nslaves = 0; endif if nslaves > 0 - global NSLAVES PARALLEL NEWORLD NSLAVES TAG; + global NSLAVES PARALLEL NEWORLD NSLAVES TAG; nn = floor(n/(NSLAVES + 1)); # number of obsns per slave # The command that the slave nodes will execute - cmd=['contrib = mle_obj_nodes(theta, data, model, modelargs, nn); ',... - 'MPI_Send(contrib,0,TAG,NEWORLD);']; + cmd=['contrib = mle_obj_nodes(theta, data, model, modelargs, nn); ',... + 'MPI_Send(contrib,0,TAG,NEWORLD);']; # send items to slaves NumCmds_Send({"theta", "nn", "cmd"}, {theta, nn, cmd}); # evaluate last block on master while slaves are busy - obj_value = mle_obj_nodes(theta, data, model, modelargs, nn); + obj_value = mle_obj_nodes(theta, data, model, modelargs, nn); # collect slaves' results contrib = 0.0; # must be initialized to use MPI_Recv - for i = 1:NSLAVES - MPI_Recv(contrib,i,TAG,NEWORLD); + for i = 1:NSLAVES + MPI_Recv(contrib,i,TAG,NEWORLD); obj_value = obj_value + contrib; endfor # compute the average - obj_value = - obj_value / n; - score = "na"; # fix this later to allow analytic score in parallel + obj_value = - obj_value / n; + score = "na"; # fix this later to allow analytic score in parallel - else # serial version - [contribs, score] = feval(model, theta, data, modelargs); + else # serial version + [contribs, score] = feval(model, theta, data, modelargs); obj_value = - mean(contribs); - if isnumeric(score) score = - mean(score)'; endif # model passes "na" when score not available - endif - + if isnumeric(score) score = - mean(score)'; endif # model passes "na" when score not available + endif # let's bullet-proof this in case the model goes nuts - if (((abs(obj_value) == Inf)) || (isnan(obj_value))) - obj_value = realmax/10; - endif + if (((abs(obj_value) == Inf)) || (isnan(obj_value))) + obj_value = realmax/10; + endif endfunction Index: nls_obj_nodes.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/econometrics/nls_obj_nodes.m,v retrieving revision 1.2 retrieving revision 1.3 diff -u -d -r1.2 -r1.3 --- nls_obj_nodes.m 25 May 2005 13:04:33 -0000 1.2 +++ nls_obj_nodes.m 30 May 2005 15:16:12 -0000 1.3 @@ -18,7 +18,7 @@ function contrib = nls_obj_nodes(theta, data, model, modelargs, nn) - global NEWORLD NSLAVES + global NEWORLD NSLAVES # Who am I? [info, rank] = MPI_Comm_rank(NEWORLD); @@ -31,7 +31,7 @@ endif data = data(startblock:endblock,:); - contrib = feval(model, theta, data, modelargs); + contrib = feval(model, theta, data, modelargs); contrib = sum(contrib); endfunction Index: gmm_estimate.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/econometrics/gmm_estimate.m,v retrieving revision 1.4 retrieving revision 1.5 diff -u -d -r1.4 -r1.5 --- gmm_estimate.m 25 May 2005 13:04:33 -0000 1.4 +++ gmm_estimate.m 30 May 2005 15:16:11 -0000 1.5 @@ -18,22 +18,26 @@ # gmm_estimate(theta, data, weight, moments, momentargs, control, nslaves) # # inputs: -# theta: column vector initial parameters -# data: data matrix -# weight: the GMM weight matrix -# moments: name of function computes the moments (should return nXg matrix of contributions) -# momentargs: (cell) additional inputs needed to compute moments. May be empty ("") -# names: vector of parameter names, e.g., use names = str2mat("param1", "param2"); -# title: string, describes model estimated -# unscale: (optional) cell that holds means and std. dev. of data (see scale_data) -# control: (optional) BFGS or SA controls (see bfgsmin and samin). May be empty (""). -# nslaves: (optional) number of slaves if executed in parallel (requires MPITB) +# theta: column vector initial parameters +# data: data matrix +# weight: the GMM weight matrix +# moments: name of function computes the moments +# (should return nXg matrix of contributions) +# momentargs: (cell) additional inputs needed to compute moments. +# May be empty ("") +# control: (optional) BFGS or SA controls (see bfgsmin and samin). +# May be empty (""). +# nslaves: (optional) number of slaves if executed in parallel +# (requires MPITB) # # outputs: # theta: GMM estimate of parameters # obj_value: the value of the gmm obj. function -# convergence: return code from bfgsmin (1 means success, see bfgsmin for details) -# iters: number of BFGS iteration used +# convergence: return code from bfgsmin +# (1 means success, see bfgsmin for details) +# iters: number of BFGS iteration used +# +# please type "gmm_example" while in octave to see an example # call the minimizing routine function [theta, obj_value, convergence, iters] = gmm_estimate(theta, data, weight, moments, momentargs, control, nslaves) @@ -62,9 +66,9 @@ if strcmp(method, "bfgs") - [theta, obj_value, convergence, iters] = bfgsmin("gmm_obj", {theta, data, weight, moments, momentargs}, control); + [theta, obj_value, convergence, iters] = bfgsmin("gmm_obj", {theta, data, weight, moments, momentargs}, control); elseif strcmp(method, "sa") - [theta, obj_value, convergence] = samin("gmm_obj", {theta, data, weight, moments, momentargs}, control); + [theta, obj_value, convergence] = samin("gmm_obj", {theta, data, weight, moments, momentargs}, control); endif if nslaves > 0 LAM_Finalize; endif # clean up Index: sum_moments_nodes.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/econometrics/sum_moments_nodes.m,v retrieving revision 1.2 retrieving revision 1.3 diff -u -d -r1.2 -r1.3 --- sum_moments_nodes.m 25 May 2005 13:04:33 -0000 1.2 +++ sum_moments_nodes.m 30 May 2005 15:16:12 -0000 1.3 @@ -18,7 +18,7 @@ function contrib = sum_moments_nodes(theta, data, moments, momentargs, nn) - global NEWORLD NSLAVES + global NEWORLD NSLAVES # Who am I? [info, rank] = MPI_Comm_rank(NEWORLD); @@ -32,7 +32,7 @@ endif data = data(startblock:endblock,:); - contrib = feval(moments, theta, data, momentargs); - contrib = sum(contrib); + contrib = feval(moments, theta, data, momentargs); + contrib = sum(contrib); endfunction Index: mle_variance.m =================================================================== RCS file: /cvsroot/octave/octave-forge/main/econometrics/mle_variance.m,v retrieving revision 1.3 retrieving revision 1.4 diff -u -d -r1.3 -r1.4 --- mle_variance.m 25 May 2005 13:04:33 -0000 1.3 +++ mle_variance.m 30 May 2005 15:16:12 -0000 1.4 @@ -21,10 +21,10 @@ # 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; + 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 |