Mercurial > forge
changeset 2360:ff994a6ff804 octave-forge
Altered the directory structure to match the package system
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/COPYING Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,339 @@ + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc. + 675 Mass Ave, Cambridge, MA 02139, USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +License is intended to guarantee your freedom to share and change free +software--to make sure the software is free for all its users. This +General Public License applies to most of the Free Software +Foundation's software and to any other program whose authors commit to +using it. (Some other Free Software Foundation software is covered by +the GNU Library General Public License instead.) You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +this service if you wish), that you receive source code or can get it +if you want it, that you can change the software or use pieces of it +in new free programs; and that you know you can do these things. + + To protect your rights, we need to make restrictions that forbid +anyone to deny you these rights or to ask you to surrender the rights. +These restrictions translate to certain responsibilities for you if you +distribute copies of the software, or if you modify it. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must give the recipients all the rights that +you have. You must make sure that they, too, receive or can get the +source code. And you must show them these terms so they know their +rights. + + We protect your rights with two steps: (1) copyright the software, and +(2) offer you this license which gives you legal permission to copy, +distribute and/or modify the software. + + Also, for each author's protection and ours, we want to make certain +that everyone understands that there is no warranty for this free +software. If the software is modified by someone else and passed on, we +want its recipients to know that what they have is not the original, so +that any problems introduced by others will not reflect on the original +authors' reputations. + + Finally, any free program is threatened constantly by software +patents. We wish to avoid the danger that redistributors of a free +program will individually obtain patent licenses, in effect making the +program proprietary. To prevent this, we have made it clear that any +patent must be licensed for everyone's free use or not licensed at all. + + The precise terms and conditions for copying, distribution and +modification follow. + + GNU GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License applies to any program or other work which contains +a notice placed by the copyright holder saying it may be distributed +under the terms of this General Public License. The "Program", below, +refers to any such program or work, and a "work based on the Program" +means either the Program or any derivative work under copyright law: +that is to say, a work containing the Program or a portion of it, +either verbatim or with modifications and/or translated into another +language. (Hereinafter, translation is included without limitation in +the term "modification".) Each licensee is addressed as "you". + +Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running the Program is not restricted, and the output from the Program +is covered only if its contents constitute a work based on the +Program (independent of having been made by running the Program). +Whether that is true depends on what the Program does. + + 1. You may copy and distribute verbatim copies of the Program's +source code as you receive it, in any medium, provided that you +conspicuously and appropriately publish on each copy an appropriate +copyright notice and disclaimer of warranty; keep intact all the +notices that refer to this License and to the absence of any warranty; +and give any other recipients of the Program a copy of this License +along with the Program. + +You may charge a fee for the physical act of transferring a copy, and +you may at your option offer warranty protection in exchange for a fee. + + 2. You may modify your copy or copies of the Program or any portion +of it, thus forming a work based on the Program, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) You must cause the modified files to carry prominent notices + stating that you changed the files and the date of any change. + + b) You must cause any work that you distribute or publish, that in + whole or in part contains or is derived from the Program or any + part thereof, to be licensed as a whole at no charge to all third + parties under the terms of this License. + + c) If the modified program normally reads commands interactively + when run, you must cause it, when started running for such + interactive use in the most ordinary way, to print or display an + announcement including an appropriate copyright notice and a + notice that there is no warranty (or else, saying that you provide + a warranty) and that users may redistribute the program under + these conditions, and telling the user how to view a copy of this + License. (Exception: if the Program itself is interactive but + does not normally print such an announcement, your work based on + the Program is not required to print an announcement.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Program, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Program, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Program. + +In addition, mere aggregation of another work not based on the Program +with the Program (or with a work based on the Program) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may copy and distribute the Program (or a work based on it, +under Section 2) in object code or executable form under the terms of +Sections 1 and 2 above provided that you also do one of the following: + + a) Accompany it with the complete corresponding machine-readable + source code, which must be distributed under the terms of Sections + 1 and 2 above on a medium customarily used for software interchange; or, + + b) Accompany it with a written offer, valid for at least three + years, to give any third party, for a charge no more than your + cost of physically performing source distribution, a complete + machine-readable copy of the corresponding source code, to be + distributed under the terms of Sections 1 and 2 above on a medium + customarily used for software interchange; or, + + c) Accompany it with the information you received as to the offer + to distribute corresponding source code. (This alternative is + allowed only for noncommercial distribution and only if you + received the program in object code or executable form with such + an offer, in accord with Subsection b above.) + +The source code for a work means the preferred form of the work for +making modifications to it. For an executable work, complete source +code means all the source code for all modules it contains, plus any +associated interface definition files, plus the scripts used to +control compilation and installation of the executable. However, as a +special exception, the source code distributed need not include +anything that is normally distributed (in either source or binary +form) with the major components (compiler, kernel, and so on) of the +operating system on which the executable runs, unless that component +itself accompanies the executable. + +If distribution of executable or object code is made by offering +access to copy from a designated place, then offering equivalent +access to copy the source code from the same place counts as +distribution of the source code, even though third parties are not +compelled to copy the source along with the object code. + + 4. You may not copy, modify, sublicense, or distribute the Program +except as expressly provided under this License. Any attempt +otherwise to copy, modify, sublicense or distribute the Program is +void, and will automatically terminate your rights under this License. +However, parties who have received copies, or rights, from you under +this License will not have their licenses terminated so long as such +parties remain in full compliance. + + 5. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Program or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Program (or any work based on the +Program), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Program or works based on it. + + 6. Each time you redistribute the Program (or any work based on the +Program), the recipient automatically receives a license from the +original licensor to copy, distribute or modify the Program subject to +these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties to +this License. + + 7. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Program at all. For example, if a patent +license would not permit royalty-free redistribution of the Program by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Program. + +If any portion of this section is held invalid or unenforceable under +any particular circumstance, the balance of the section is intended to +apply and the section as a whole is intended to apply in other +circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system, which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 8. If the distribution and/or use of the Program is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Program under this License +may add an explicit geographical distribution limitation excluding +those countries, so that distribution is permitted only in or among +countries not thus excluded. In such case, this License incorporates +the limitation as if written in the body of this License. + + 9. The Free Software Foundation may publish revised and/or new versions +of the General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + +Each version is given a distinguishing version number. If the Program +specifies a version number of this License which applies to it and "any +later version", you have the option of following the terms and conditions +either of that version or of any later version published by the Free +Software Foundation. If the Program does not specify a version number of +this License, you may choose any version ever published by the Free Software +Foundation. + + 10. If you wish to incorporate parts of the Program into other free +programs whose distribution conditions are different, write to the author +to ask for permission. For software which is copyrighted by the Free +Software Foundation, write to the Free Software Foundation; we sometimes +make exceptions for this. Our decision will be guided by the two goals +of preserving the free status of all derivatives of our free software and +of promoting the sharing and reuse of software generally. + + NO WARRANTY + + 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY +FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN +OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES +PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED +OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS +TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE +PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, +REPAIR OR CORRECTION. + + 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR +REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, +INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING +OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED +TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY +YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER +PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE +POSSIBILITY OF SUCH DAMAGES. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +convey the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + <one line to give the program's name and a brief idea of what it does.> + Copyright (C) 19yy <name of author> + + 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., 675 Mass Ave, Cambridge, MA 02139, USA. + +Also add information on how to contact you by electronic and paper mail. + +If the program is interactive, make it output a short notice like this +when it starts in an interactive mode: + + Gnomovision version 69, Copyright (C) 19yy name of author + Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, the commands you use may +be called something other than `show w' and `show c'; they could even be +mouse-clicks or menu items--whatever suits your program. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the program, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the program + `Gnomovision' (which makes passes at compilers) written by James Hacker. + + <signature of Ty Coon>, 1 April 1989 + Ty Coon, President of Vice + +This General Public License does not permit incorporating your program into +proprietary programs. If your program is a subroutine library, you may +consider it more useful to permit linking proprietary applications with the +library. If this is what you want to do, use the GNU Library General +Public License instead of this License.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/DESCRIPTION Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,10 @@ +Name: Econometrics +Version: 1.0.0 +Date: 2006-08-05 +Author: Michael Creel <michael.creel@uab.es> +Maintainer: Michael Creel <michael.creel@uab.es> +Title: Econometrics. +Description: Add a description to this package! +Depends: octave (>= 2.9.7) +License: GPL version 2 or later +Url: http://octave.sf.net
--- a/main/econometrics/README Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,14 +0,0 @@ -The main functions for end users are: - -mle_estimate.m Performs estimation by MLE -mle_results.m Performs estimatiion by MLE and shows results -mle_example.m Example - - -gmm_estimate.m Performs estimation by GMM -gmm_results.m Performs estimatiion by GMM and shows results -gmm_example.m Example - - -The other functions are mostly for internal use, or for use by -advanced users.
--- a/main/econometrics/average_moments.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,55 +0,0 @@ -# Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es> -# -# 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 - -# for internal use by gmm_estimate - - -# average moments (separate function so it can be differentiated) -function m = average_moments(theta, data, moments, momentargs) - - global NSLAVES PARALLEL NEWORLD NSLAVES TAG; - - n = rows(data); - - 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);']; - - # 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); - - # collect slaves' results - contrib = zeros(1,columns(m)); - for i = 1:NSLAVES - MPI_Recv(contrib,i,TAG,NEWORLD); - m = m + contrib; - endfor - - m = m'; # we want a column vector, please - m = m/n; # average please, not sum - - else # serial version - m = feval(moments, theta, data, momentargs); - m = mean(m)'; # returns Gx1 moment vector - endif - -endfunction
--- a/main/econometrics/delta_method.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,25 +0,0 @@ -# Copyright (C) 2003,2004 Michael Creel <michael.creel@uab.es> -# -# 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, var_transf] = delta_method(func, theta, otherargs, vartheta) - theta_transf = feval(func, theta, otherargs); - D = numgradient(func, {theta, otherargs}); - - var_transf = D * vartheta * D'; - -endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/doc/README Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,14 @@ +The main functions for end users are: + +mle_estimate.m Performs estimation by MLE +mle_results.m Performs estimatiion by MLE and shows results +mle_example.m Example + + +gmm_estimate.m Performs estimation by GMM +gmm_results.m Performs estimatiion by GMM and shows results +gmm_example.m Example + + +The other functions are mostly for internal use, or for use by +advanced users.
--- a/main/econometrics/gmm_estimate.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,76 +0,0 @@ -# Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es> -# -# 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 - -# usage: [theta, obj_value, convergence, iters] = -# 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 ("") -# 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 -# -# 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) - - if nargin < 5 error("gmm_estimate: 5 arguments required"); endif - - if nargin < 6 control = {Inf,0,1,1}; endif # default controls - if !iscell(control) control = {Inf,0,1,1}; endif # default controls if receive placeholder - if nargin < 7 nslaves = 0; endif - - if nslaves > 0 - global NSLAVES PARALLEL NEWORLD NSLAVES TAG; - LAM_Init(nslaves); - # Send the data to all nodes - NumCmds_Send({"data", "weight", "moments", "momentargs"}, {data, weight, moments, momentargs}); - endif - - # bfgs or sa? - if (size(control,1)*size(control,2) == 0) # use default bfgs if no control - control = {Inf,0,1,1}; - method = "bfgs"; - elseif (size(control,1)*size(control,2) < 11) - method = "bfgs"; - else method = "sa"; - endif - - - if strcmp(method, "bfgs") - [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); - endif - - if nslaves > 0 LAM_Finalize; endif # clean up - -endfunction
--- a/main/econometrics/gmm_example.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,65 +0,0 @@ -# Copyright (C) 2003,2004, 2005 Michael Creel <michael.creel@uab.es> -# -# 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 - -# GMM example file, shows initial consistent estimator, -# estimation of efficient weight, and second round -# efficient estimator - - -n = 1000; -k = 5; - -x = [ones(n,1) rand(n,k-1)]; -w = [x, rand(n,1)]; -theta = [-k+1; ones(k-1,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 = "poisson_moments"; -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,2,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, control); - -# 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); - - -# Example doing estimation in parallel on a cluster (requires MPITB) -# 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);
--- a/main/econometrics/gmm_obj.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,32 +0,0 @@ -# Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es> -# -# 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 - - -# The GMM objective function, for internal use by gmm_estimate -# 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
--- a/main/econometrics/gmm_results.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,107 +0,0 @@ -# Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es> -# -# 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 - -# usage: [theta, V, obj_value] = -# 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., 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 (inverse of covariance of moments) -# obj_value: the value of the GMM objective function -# -# 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) - - if nargin < 10 nslaves = 0; endif # serial by default - - if nargin < 9 - [theta, obj_value, convergence] = gmm_estimate(theta, data, weight, moments, momentargs, "", nslaves); - else - [theta, obj_value, convergence] = gmm_estimate(theta, data, weight, moments, momentargs, control, nslaves); - 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
--- a/main/econometrics/gmm_variance.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,25 +0,0 @@ -# Copyright (C) 2003,2004 Michael Creel <michael.creel@uab.es> -# -# 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 - - -# 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
--- a/main/econometrics/gmm_variance_inefficient.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,28 +0,0 @@ -# Copyright (C) 2003,2004 Michael Creel <michael.creel@uab.es> -# -# 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 - -# 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'; - - n = rows(data); - - J = D*weight*D'; - J = inv(J); - I = D*weight*omega*weight*D'; - V = (1/n)*J*I*J; -endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/average_moments.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,55 @@ +# Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es> +# +# 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 + +# for internal use by gmm_estimate + + +# average moments (separate function so it can be differentiated) +function m = average_moments(theta, data, moments, momentargs) + + global NSLAVES PARALLEL NEWORLD NSLAVES TAG; + + n = rows(data); + + 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);']; + + # 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); + + # collect slaves' results + contrib = zeros(1,columns(m)); + for i = 1:NSLAVES + MPI_Recv(contrib,i,TAG,NEWORLD); + m = m + contrib; + endfor + + m = m'; # we want a column vector, please + m = m/n; # average please, not sum + + else # serial version + m = feval(moments, theta, data, momentargs); + m = mean(m)'; # returns Gx1 moment vector + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/delta_method.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,25 @@ +# Copyright (C) 2003,2004 Michael Creel <michael.creel@uab.es> +# +# 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, var_transf] = delta_method(func, theta, otherargs, vartheta) + theta_transf = feval(func, theta, otherargs); + D = numgradient(func, {theta, otherargs}); + + var_transf = D * vartheta * D'; + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/gmm_estimate.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,76 @@ +# Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es> +# +# 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 + +# usage: [theta, obj_value, convergence, iters] = +# 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 ("") +# 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 +# +# 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) + + if nargin < 5 error("gmm_estimate: 5 arguments required"); endif + + if nargin < 6 control = {Inf,0,1,1}; endif # default controls + if !iscell(control) control = {Inf,0,1,1}; endif # default controls if receive placeholder + if nargin < 7 nslaves = 0; endif + + if nslaves > 0 + global NSLAVES PARALLEL NEWORLD NSLAVES TAG; + LAM_Init(nslaves); + # Send the data to all nodes + NumCmds_Send({"data", "weight", "moments", "momentargs"}, {data, weight, moments, momentargs}); + endif + + # bfgs or sa? + if (size(control,1)*size(control,2) == 0) # use default bfgs if no control + control = {Inf,0,1,1}; + method = "bfgs"; + elseif (size(control,1)*size(control,2) < 11) + method = "bfgs"; + else method = "sa"; + endif + + + if strcmp(method, "bfgs") + [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); + endif + + if nslaves > 0 LAM_Finalize; endif # clean up + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/gmm_example.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,65 @@ +# Copyright (C) 2003,2004, 2005 Michael Creel <michael.creel@uab.es> +# +# 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 + +# GMM example file, shows initial consistent estimator, +# estimation of efficient weight, and second round +# efficient estimator + + +n = 1000; +k = 5; + +x = [ones(n,1) rand(n,k-1)]; +w = [x, rand(n,1)]; +theta = [-k+1; ones(k-1,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 = "poisson_moments"; +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,2,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, control); + +# 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); + + +# Example doing estimation in parallel on a cluster (requires MPITB) +# 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);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/gmm_obj.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,32 @@ +# Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es> +# +# 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 + + +# The GMM objective function, for internal use by gmm_estimate +# 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/gmm_results.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,107 @@ +# Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es> +# +# 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 + +# usage: [theta, V, obj_value] = +# 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., 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 (inverse of covariance of moments) +# obj_value: the value of the GMM objective function +# +# 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) + + if nargin < 10 nslaves = 0; endif # serial by default + + if nargin < 9 + [theta, obj_value, convergence] = gmm_estimate(theta, data, weight, moments, momentargs, "", nslaves); + else + [theta, obj_value, convergence] = gmm_estimate(theta, data, weight, moments, momentargs, control, nslaves); + 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/gmm_variance.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,25 @@ +# Copyright (C) 2003,2004 Michael Creel <michael.creel@uab.es> +# +# 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 + + +# 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/gmm_variance_inefficient.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,28 @@ +# Copyright (C) 2003,2004 Michael Creel <michael.creel@uab.es> +# +# 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 + +# 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'; + + n = rows(data); + + J = D*weight*D'; + J = inv(J); + I = D*weight*omega*weight*D'; + V = (1/n)*J*I*J; +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/mle_estimate.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,74 @@ +## Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es> +## +## 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 + +# usage: +# [theta, obj_value, conv, iters] = mle_estimate(theta, data, model, modelargs, control, nslaves) +# +# inputs: +# theta: column vector of model parameters +# data: data matrix +# model: name of function that computes log-likelihood +# modelargs: (cell) additional inputs needed by model. 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: ML estimated value of parameters +# obj_value: the value of the log likelihood function at ML estimate +# conv: return code from bfgsmin (1 means success, see bfgsmin for details) +# iters: number of BFGS iteration used +# +# please see mle_example.m for examples of how to use this +function [theta, obj_value, convergence, iters] = mle_estimate(theta, data, model, modelargs, control, nslaves) + + + if nargin < 3 + error("mle_estimate: 3 arguments required"); + endif + + if nargin < 4 modelargs = {}; endif # create placeholder if not used + if !iscell(modelargs) modelargs = {}; endif # default controls if receive placeholder + if nargin < 5 control = {Inf,0,1,1}; endif # default controls and method + if !iscell(control) control = {Inf,0,1,1}; endif # default controls if receive placeholder + if nargin < 6 nslaves = 0; endif + if nslaves > 0 + global NSLAVES PARALLEL NEWORLD TAG; + LAM_Init(nslaves); + # Send the data to all nodes + NumCmds_Send({"data", "model", "modelargs"}, {data, model, modelargs}); + endif + + # bfgs or sa? + if (size(control,1)*size(control,2) == 0) # use default bfgs if no control + control = {Inf,0,1,1}; + method = "bfgs"; + elseif (size(control,1)*size(control,2) < 11) + method = "bfgs"; + else method = "sa"; + endif + + + if strcmp(method, "bfgs") + [theta, obj_value, convergence, iters] = bfgsmin("mle_obj", {theta, data, model, modelargs, nslaves}, control); + elseif strcmp(method, "sa") + [theta, obj_value, convergence] = samin("mle_obj", {theta, data, model, modelargs, nslaves}, control); + endif + + if nslaves > 0 + LAM_Finalize; + endif # cleanup + obj_value = - obj_value; # recover from minimization rather than maximization +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/mle_example.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,80 @@ +## Copyright (C) 2003,2004 Michael Creel <michael.creel@uab.es> +## +## 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 + + + +# 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"; +modelargs = {0}; # if this is zero the function gives analytic score, otherwise not +# 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 = {50,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); +modelargs = {1}; # set the switch for no score +[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) +# 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);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/mle_obj.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,65 @@ +# Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es> +# +# 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 + +## usage: [obj_value, score] = mle_obj(theta, data, model, modelargs, nslaves) +## +## Returns the average log-likelihood for a specified model +## This is for internal use by mle_estimate + + +function [obj_value, score] = mle_obj(theta, data, model, modelargs, nslaves) + + n = rows(data); + + if nargin < 5 nslaves = 0; endif + if nslaves > 0 + 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);']; + + # 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); + + # collect slaves' results + contrib = 0.0; # must be initialized to use MPI_Recv + 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 + + 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 + + # let's bullet-proof this in case the model goes nuts + if (((abs(obj_value) == Inf)) || (isnan(obj_value))) + obj_value = realmax/10; + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/mle_obj_nodes.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,34 @@ +## Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es> +## +## 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 + +function contrib = mle_obj_nodes(theta, data, model, modelargs, nn) + global NEWORLD NSLAVES + + # Who am I? + [info, rank] = MPI_Comm_rank(NEWORLD); + if rank == 0 # Do this if I'm master + startblock = NSLAVES*nn + 1; + endblock = rows(data); + else # this is for the slaves + startblock = rank*nn-nn+1; + endblock = rank*nn; + endif + + data = data(startblock:endblock,:); + contrib = feval(model, theta, data, modelargs); + contrib = sum(contrib); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/mle_results.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,93 @@ +# Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es> +# +# 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 + +# usage: [theta, V, obj_value, infocrit] = +# mle_results(theta, data, model, modelargs, names, title, unscale, control) +# +# inputs: +# theta: column vector of model parameters +# data: data matrix +# model: name of function that computes log-likelihood +# modelargs: (cell) additional inputs needed by model. 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) +# +# outputs: +# theta: ML estimated value of parameters +# obj_value: the value of the log likelihood function at ML estimate +# conv: return code from bfgsmin (1 means success, see bfgsmin for details) +# iters: number of BFGS iteration used + + +## +## Please see mle_example for information on how to use this + +# 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 < 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); + + # unscale results if argument has been passed + # this puts coefficients into scale corresponding to the original modelargs + if (nargin > 6) + 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/mle_variance.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,30 @@ +## Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es> +## +## 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 + +## usage: [V,scorecontribs,J_inv] = +## mle_variance(theta, data, model, modelargs) +## +## This is for internal use by mle_results + +# 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/nls_estimate.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,71 @@ +## Copyright (C) 2005 Michael Creel <michael.creel@uab.es> +## +## 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 + +# usage: +# [theta, obj_value, conv, iters] = nls_estimate(theta, data, model, modelargs, control, nslaves) +# +# inputs: +# theta: column vector of model parameters +# data: data matrix +# model: name of function that computes the vector of sums of squared errors +# modelargs: (cell) additional inputs needed by model. 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: NLS estimated value of parameters +# obj_value: the value of the sum of squared errors at NLS estimate +# conv: return code from bfgsmin (1 means success, see bfgsmin for details) +# iters: number of BFGS iteration used +# +# please see nls_example.m for examples of how to use this +function [theta, obj_value, convergence, iters] = nls_estimate(theta, data, model, modelargs, control, nslaves) + + if nargin < 3 + error("nls_estimate: 3 arguments required"); + endif + + if nargin < 4 modelargs = {}; endif # create placeholder if not used + if !iscell(modelargs) modelargs = {}; endif # default controls if receive placeholder + if nargin < 5 control = {Inf,0,1,1}; endif # default controls and method + if !iscell(control) control = {Inf,0,1,1}; endif # default controls if receive placeholder + if nargin < 6 nslaves = 0; endif + + if nslaves > 0 + global NSLAVES PARALLEL NEWORLD TAG + LAM_Init(nslaves); + # Send the data to all nodes + NumCmds_Send({"data", "model", "modelargs"}, {data, model, modelargs}); + endif + + # bfgs or sa? + if (size(control,1)*size(control,2) == 0) # use default bfgs if no control + control = {Inf,0,1,1}; + method = "bfgs"; + elseif (size(control,1)*size(control,2) < 11) + method = "bfgs"; + else method = "sa"; + endif + + if strcmp(method, "bfgs") + [theta, obj_value, convergence, iters] = bfgsmin("nls_obj", {theta, data, model, modelargs, nslaves}, control); + elseif strcmp(method, "sa") + [theta, obj_value, convergence] = samin("nls_obj", {theta, data, model, modelargs, nslaves}, control); + endif + + # cleanup + if nslaves > 0 LAM_Finalize; endif +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/nls_example.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,60 @@ +## Copyright (C) 2006 Michael Creel <michael.creel@uab.es> +## +## 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 NLS + +# Generate data +n = 100; # how many observations? + +# the explanatory variables: note that they have unequal scales +x = [ones(n,1) rand(n,2)]; +theta = 1:3; # true coefficients are 1,2,3 +theta = theta'; + +lambda = exp(x*theta); +y = randp(lambda); # generate the dependent variable + +# example objective function for nls +function [obj_contrib, score] = nls_example_obj(theta, data, otherargs) + y = data(:,1); + x = data(:,2:columns(data)); + lambda = exp(x*theta); + errors = y - lambda; + obj_contrib = errors .* errors; + score = "na"; +endfunction + + +#################################### +# define arguments for nls_estimate # +#################################### + +# starting values +theta = zeros(3,1); + +# data +data = [y, x]; + +# name of model to estimate +model = "nls_example_obj"; +modelargs = {}; # none required for this obj fn. + + +# controls for bfgsmin +control = {50,1,1,1}; + +printf("\n\NLS estimation example\n"); +[theta, obj_value, convergence] = nls_estimate(theta, data, model, modelargs, control);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/nls_obj.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,61 @@ +# Copyright (C) 2005 Michael Creel <michael.creel@uab.es> +# +# 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 + +# usage: [obj_value, score] = nls_obj(theta, data, model, modelargs, nslaves) +# +# Returns the average sum of squared errors for a specified model +# This is for internal use by nls_estimate + + +function [obj_value, score] = nls_obj(theta, data, model, modelargs, nslaves) + + n = rows(data); + + 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);']; + + # 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); + + # collect slaves' results + contrib = 0.0; # must be initialized to use MPI_Recv + 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 + + 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 + + # let's bullet-proof this in case the model goes nuts + if (((abs(obj_value) == Inf)) || (isnan(obj_value))) obj_value = realmax; endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/nls_obj_nodes.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,37 @@ +# Copyright (C) 2005 Michael Creel <michael.creel@uab.es> +# +# 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 for internal use by nls_estimate + + +function contrib = nls_obj_nodes(theta, data, model, modelargs, nn) + global NEWORLD NSLAVES + # Who am I? + [info, rank] = MPI_Comm_rank(NEWORLD); + + if rank == 0 # Do this if I'm master + startblock = NSLAVES*nn + 1; + endblock = rows(data); + else # this is for the slaves + startblock = rank*nn-nn+1; + endblock = rank*nn; + endif + + data = data(startblock:endblock,:); + contrib = feval(model, theta, data, modelargs); + contrib = sum(contrib); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/parameterize.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,29 @@ +# Copyright (C) 2003,2004 Michael Creel <michael.creel@uab.es> +# +# 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 +# + +# usage: theta = parameterize(theta, otherargs) +# +# 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/poisson.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,26 @@ +# Copyright (C) 2005 Michael Creel <michael.creel@uab.es> +# +# 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 likelihood function (Poisson for count data) with score + +function [log_density, score] = poisson(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); + if (otherargs{1} == 1) score = "na"; endif # provide analytic score or not? +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/poisson_moments.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,27 @@ +# Copyright (C) 2005 Michael Creel <michael.creel@uab.es> +# +# 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 + +# the form a user-written moment function should take + +function m = poisson_moments(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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/prettyprint.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,48 @@ +# Copyright (C) 2003,2004 Michael Creel <michael.creel@uab.es> +# +# 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 ischar(rlabels(i,:)) + printf(rlabels(i,:)); + else + printf("%i", rlabels(i,:)); + endif + printf(" %10.3f", mat(i,:)); + printf("\n"); + endfor +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/prettyprint_c.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,38 @@ +# Copyright (C) 2003,2004 Michael Creel <michael.creel@uab.es> +# +# 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/scale_data.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,47 @@ +# Copyright (C) 2003,2004 Michael Creel <michael.creel@uab.es> +# +# 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/sum_moments_nodes.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,38 @@ +# Copyright (C) 2005 Michael Creel <michael.creel@uab.es> +# +# 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 + +# for internal use by gmm_estimate + +function contrib = sum_moments_nodes(theta, data, moments, momentargs, nn) + + global NEWORLD NSLAVES + + # Who am I? + [info, rank] = MPI_Comm_rank(NEWORLD); + + if rank == 0 # Do this if I'm master + startblock = NSLAVES*nn + 1; + endblock = rows(data); + else # this is for the slaves + startblock = rank*nn-nn+1; + endblock = rank*nn; + endif + + data = data(startblock:endblock,:); + contrib = feval(moments, theta, data, momentargs); + contrib = sum(contrib); + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/econometrics/inst/unscale_parameters.m Sun Aug 20 12:18:19 2006 +0000 @@ -0,0 +1,40 @@ +# Copyright (C) 2003,2004 Michael Creel <michael.creel@uab.es> +# +# 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
--- a/main/econometrics/mle_estimate.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,74 +0,0 @@ -## Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es> -## -## 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 - -# usage: -# [theta, obj_value, conv, iters] = mle_estimate(theta, data, model, modelargs, control, nslaves) -# -# inputs: -# theta: column vector of model parameters -# data: data matrix -# model: name of function that computes log-likelihood -# modelargs: (cell) additional inputs needed by model. 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: ML estimated value of parameters -# obj_value: the value of the log likelihood function at ML estimate -# conv: return code from bfgsmin (1 means success, see bfgsmin for details) -# iters: number of BFGS iteration used -# -# please see mle_example.m for examples of how to use this -function [theta, obj_value, convergence, iters] = mle_estimate(theta, data, model, modelargs, control, nslaves) - - - if nargin < 3 - error("mle_estimate: 3 arguments required"); - endif - - if nargin < 4 modelargs = {}; endif # create placeholder if not used - if !iscell(modelargs) modelargs = {}; endif # default controls if receive placeholder - if nargin < 5 control = {Inf,0,1,1}; endif # default controls and method - if !iscell(control) control = {Inf,0,1,1}; endif # default controls if receive placeholder - if nargin < 6 nslaves = 0; endif - if nslaves > 0 - global NSLAVES PARALLEL NEWORLD TAG; - LAM_Init(nslaves); - # Send the data to all nodes - NumCmds_Send({"data", "model", "modelargs"}, {data, model, modelargs}); - endif - - # bfgs or sa? - if (size(control,1)*size(control,2) == 0) # use default bfgs if no control - control = {Inf,0,1,1}; - method = "bfgs"; - elseif (size(control,1)*size(control,2) < 11) - method = "bfgs"; - else method = "sa"; - endif - - - if strcmp(method, "bfgs") - [theta, obj_value, convergence, iters] = bfgsmin("mle_obj", {theta, data, model, modelargs, nslaves}, control); - elseif strcmp(method, "sa") - [theta, obj_value, convergence] = samin("mle_obj", {theta, data, model, modelargs, nslaves}, control); - endif - - if nslaves > 0 - LAM_Finalize; - endif # cleanup - obj_value = - obj_value; # recover from minimization rather than maximization -endfunction
--- a/main/econometrics/mle_example.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,80 +0,0 @@ -## Copyright (C) 2003,2004 Michael Creel <michael.creel@uab.es> -## -## 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 - - - -# 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"; -modelargs = {0}; # if this is zero the function gives analytic score, otherwise not -# 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 = {50,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); -modelargs = {1}; # set the switch for no score -[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) -# 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);
--- a/main/econometrics/mle_obj.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,65 +0,0 @@ -# Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es> -# -# 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 - -## usage: [obj_value, score] = mle_obj(theta, data, model, modelargs, nslaves) -## -## Returns the average log-likelihood for a specified model -## This is for internal use by mle_estimate - - -function [obj_value, score] = mle_obj(theta, data, model, modelargs, nslaves) - - n = rows(data); - - if nargin < 5 nslaves = 0; endif - if nslaves > 0 - 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);']; - - # 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); - - # collect slaves' results - contrib = 0.0; # must be initialized to use MPI_Recv - 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 - - 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 - - # let's bullet-proof this in case the model goes nuts - if (((abs(obj_value) == Inf)) || (isnan(obj_value))) - obj_value = realmax/10; - endif - -endfunction
--- a/main/econometrics/mle_obj_nodes.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,34 +0,0 @@ -## Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es> -## -## 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 - -function contrib = mle_obj_nodes(theta, data, model, modelargs, nn) - global NEWORLD NSLAVES - - # Who am I? - [info, rank] = MPI_Comm_rank(NEWORLD); - if rank == 0 # Do this if I'm master - startblock = NSLAVES*nn + 1; - endblock = rows(data); - else # this is for the slaves - startblock = rank*nn-nn+1; - endblock = rank*nn; - endif - - data = data(startblock:endblock,:); - contrib = feval(model, theta, data, modelargs); - contrib = sum(contrib); - -endfunction
--- a/main/econometrics/mle_results.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,93 +0,0 @@ -# Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es> -# -# 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 - -# usage: [theta, V, obj_value, infocrit] = -# mle_results(theta, data, model, modelargs, names, title, unscale, control) -# -# inputs: -# theta: column vector of model parameters -# data: data matrix -# model: name of function that computes log-likelihood -# modelargs: (cell) additional inputs needed by model. 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) -# -# outputs: -# theta: ML estimated value of parameters -# obj_value: the value of the log likelihood function at ML estimate -# conv: return code from bfgsmin (1 means success, see bfgsmin for details) -# iters: number of BFGS iteration used - - -## -## Please see mle_example for information on how to use this - -# 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 < 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); - - # unscale results if argument has been passed - # this puts coefficients into scale corresponding to the original modelargs - if (nargin > 6) - 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
--- a/main/econometrics/mle_variance.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,30 +0,0 @@ -## Copyright (C) 2003,2004,2005 Michael Creel <michael.creel@uab.es> -## -## 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 - -## usage: [V,scorecontribs,J_inv] = -## mle_variance(theta, data, model, modelargs) -## -## This is for internal use by mle_results - -# 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
--- a/main/econometrics/nls_estimate.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,71 +0,0 @@ -## Copyright (C) 2005 Michael Creel <michael.creel@uab.es> -## -## 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 - -# usage: -# [theta, obj_value, conv, iters] = nls_estimate(theta, data, model, modelargs, control, nslaves) -# -# inputs: -# theta: column vector of model parameters -# data: data matrix -# model: name of function that computes the vector of sums of squared errors -# modelargs: (cell) additional inputs needed by model. 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: NLS estimated value of parameters -# obj_value: the value of the sum of squared errors at NLS estimate -# conv: return code from bfgsmin (1 means success, see bfgsmin for details) -# iters: number of BFGS iteration used -# -# please see nls_example.m for examples of how to use this -function [theta, obj_value, convergence, iters] = nls_estimate(theta, data, model, modelargs, control, nslaves) - - if nargin < 3 - error("nls_estimate: 3 arguments required"); - endif - - if nargin < 4 modelargs = {}; endif # create placeholder if not used - if !iscell(modelargs) modelargs = {}; endif # default controls if receive placeholder - if nargin < 5 control = {Inf,0,1,1}; endif # default controls and method - if !iscell(control) control = {Inf,0,1,1}; endif # default controls if receive placeholder - if nargin < 6 nslaves = 0; endif - - if nslaves > 0 - global NSLAVES PARALLEL NEWORLD TAG - LAM_Init(nslaves); - # Send the data to all nodes - NumCmds_Send({"data", "model", "modelargs"}, {data, model, modelargs}); - endif - - # bfgs or sa? - if (size(control,1)*size(control,2) == 0) # use default bfgs if no control - control = {Inf,0,1,1}; - method = "bfgs"; - elseif (size(control,1)*size(control,2) < 11) - method = "bfgs"; - else method = "sa"; - endif - - if strcmp(method, "bfgs") - [theta, obj_value, convergence, iters] = bfgsmin("nls_obj", {theta, data, model, modelargs, nslaves}, control); - elseif strcmp(method, "sa") - [theta, obj_value, convergence] = samin("nls_obj", {theta, data, model, modelargs, nslaves}, control); - endif - - # cleanup - if nslaves > 0 LAM_Finalize; endif -endfunction
--- a/main/econometrics/nls_example.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,60 +0,0 @@ -## Copyright (C) 2006 Michael Creel <michael.creel@uab.es> -## -## 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 NLS - -# Generate data -n = 100; # how many observations? - -# the explanatory variables: note that they have unequal scales -x = [ones(n,1) rand(n,2)]; -theta = 1:3; # true coefficients are 1,2,3 -theta = theta'; - -lambda = exp(x*theta); -y = randp(lambda); # generate the dependent variable - -# example objective function for nls -function [obj_contrib, score] = nls_example_obj(theta, data, otherargs) - y = data(:,1); - x = data(:,2:columns(data)); - lambda = exp(x*theta); - errors = y - lambda; - obj_contrib = errors .* errors; - score = "na"; -endfunction - - -#################################### -# define arguments for nls_estimate # -#################################### - -# starting values -theta = zeros(3,1); - -# data -data = [y, x]; - -# name of model to estimate -model = "nls_example_obj"; -modelargs = {}; # none required for this obj fn. - - -# controls for bfgsmin -control = {50,1,1,1}; - -printf("\n\NLS estimation example\n"); -[theta, obj_value, convergence] = nls_estimate(theta, data, model, modelargs, control);
--- a/main/econometrics/nls_obj.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,61 +0,0 @@ -# Copyright (C) 2005 Michael Creel <michael.creel@uab.es> -# -# 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 - -# usage: [obj_value, score] = nls_obj(theta, data, model, modelargs, nslaves) -# -# Returns the average sum of squared errors for a specified model -# This is for internal use by nls_estimate - - -function [obj_value, score] = nls_obj(theta, data, model, modelargs, nslaves) - - n = rows(data); - - 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);']; - - # 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); - - # collect slaves' results - contrib = 0.0; # must be initialized to use MPI_Recv - 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 - - 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 - - # let's bullet-proof this in case the model goes nuts - if (((abs(obj_value) == Inf)) || (isnan(obj_value))) obj_value = realmax; endif - -endfunction
--- a/main/econometrics/nls_obj_nodes.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,37 +0,0 @@ -# Copyright (C) 2005 Michael Creel <michael.creel@uab.es> -# -# 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 for internal use by nls_estimate - - -function contrib = nls_obj_nodes(theta, data, model, modelargs, nn) - global NEWORLD NSLAVES - # Who am I? - [info, rank] = MPI_Comm_rank(NEWORLD); - - if rank == 0 # Do this if I'm master - startblock = NSLAVES*nn + 1; - endblock = rows(data); - else # this is for the slaves - startblock = rank*nn-nn+1; - endblock = rank*nn; - endif - - data = data(startblock:endblock,:); - contrib = feval(model, theta, data, modelargs); - contrib = sum(contrib); - -endfunction
--- a/main/econometrics/parameterize.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,29 +0,0 @@ -# Copyright (C) 2003,2004 Michael Creel <michael.creel@uab.es> -# -# 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 -# - -# usage: theta = parameterize(theta, otherargs) -# -# 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
--- a/main/econometrics/poisson.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,26 +0,0 @@ -# Copyright (C) 2005 Michael Creel <michael.creel@uab.es> -# -# 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 likelihood function (Poisson for count data) with score - -function [log_density, score] = poisson(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); - if (otherargs{1} == 1) score = "na"; endif # provide analytic score or not? -endfunction
--- a/main/econometrics/poisson_moments.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,27 +0,0 @@ -# Copyright (C) 2005 Michael Creel <michael.creel@uab.es> -# -# 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 - -# the form a user-written moment function should take - -function m = poisson_moments(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
--- a/main/econometrics/prettyprint.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,48 +0,0 @@ -# Copyright (C) 2003,2004 Michael Creel <michael.creel@uab.es> -# -# 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 ischar(rlabels(i,:)) - printf(rlabels(i,:)); - else - printf("%i", rlabels(i,:)); - endif - printf(" %10.3f", mat(i,:)); - printf("\n"); - endfor -endfunction
--- a/main/econometrics/prettyprint_c.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,38 +0,0 @@ -# Copyright (C) 2003,2004 Michael Creel <michael.creel@uab.es> -# -# 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
--- a/main/econometrics/scale_data.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,47 +0,0 @@ -# Copyright (C) 2003,2004 Michael Creel <michael.creel@uab.es> -# -# 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
--- a/main/econometrics/sum_moments_nodes.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,38 +0,0 @@ -# Copyright (C) 2005 Michael Creel <michael.creel@uab.es> -# -# 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 - -# for internal use by gmm_estimate - -function contrib = sum_moments_nodes(theta, data, moments, momentargs, nn) - - global NEWORLD NSLAVES - - # Who am I? - [info, rank] = MPI_Comm_rank(NEWORLD); - - if rank == 0 # Do this if I'm master - startblock = NSLAVES*nn + 1; - endblock = rows(data); - else # this is for the slaves - startblock = rank*nn-nn+1; - endblock = rank*nn; - endif - - data = data(startblock:endblock,:); - contrib = feval(moments, theta, data, momentargs); - contrib = sum(contrib); - -endfunction
--- a/main/econometrics/unscale_parameters.m Sun Aug 20 12:12:23 2006 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,40 +0,0 @@ -# Copyright (C) 2003,2004 Michael Creel <michael.creel@uab.es> -# -# 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