Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: st: RE: Maximize a function that contains an integral


From   Bob Hammond <robert_hammond@ncsu.edu>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: RE: Maximize a function that contains an integral
Date   Wed, 04 Mar 2009 14:09:58 -0500

All,

Thanks to Al for directing me to his presentation on numerical integration methods. I think that I figured out how to manipulate the bounds of integration in the Mata function nwspgr() based on Al's discussion on the Stata command integ. nwspgr() (available in a zip file at http://sparse-grids.de/#Stata) creates a matrix (call it x) and a column vector (call it w) that can be used to approximate an integral defined as a Mata function. I can manipulate the bounds of integration by manipulating the x matrix prior to multiplying it by the weighting matrix. I've pasted a sample program below with two simple integrals (that have analytical solutions) and run them through some simulations (not shown). My integral is of one dimension so I have not played much with the nwspgr()'s ability to handle integrals of higher dimensions. With one dimension, for all lower and upper bounds and a variety of integrals, nwspgr() gives approximation errors no more than 10^-7. I've also tried integrals without analytical solutions and compared to Mathematica with similarly precise approximations.

This may have been clear to those who know more about numerical integration than I but it was new to me. Take care,

Bob

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
cap mata: mata drop g()
cap mata: mata drop h()
mata:
   real colvector g(real matrix x){
       return( x:^4 )
   }
   real colvector h(real matrix x){
      return( normalden(x) )
} a = -1 // lower bound of integration
   b = 2        // upper bound of integration
   D = 1       // dimensions of integration
   k = 25      // accuracy level of approximation
   nw = nwspgr("KPU", D, k)
   x = nw[.,1..D]:*(b-a) :+ a
   w = nw[.,D+1]
   trueval = (1/5)*b^5 - (1/5)*a^5
   approx  = g(x)'*w*(b-a)
   error   = sqrt((approx - trueval):^2):/trueval
printf("True Value = %10.0g, Approximation = %10.0g, Approximation Error = %10.0g\n",trueval,approx,error)
   trueval = normal(b) - normal(a)
   approx  = h(x)'*w*(b-a)
   error   = sqrt((approx - trueval):^2):/trueval
printf("True Value = %10.0g, Approximation = %10.0g, Approximation Error = %10.0g\n",trueval,approx,error)
end
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

Feiveson, Alan H. (JSC-SK311) wrote:
Bob - Try looking at http://www.stata.com/meeting/3nasug/abstracts.html
Al Feiveson

-----Original Message-----
From: owner-statalist@hsphsun2.harvard.edu [mailto:owner-statalist@hsphsun2.harvard.edu] On Behalf Of Bob Hammond
Sent: Tuesday, March 03, 2009 4:24 PM
To: statalist@hsphsun2.harvard.edu
Subject: st: Maximize a function that contains an integral

All,

I want to find the maximum of a function that contains an integral without an analytical solution. The difficulty is that the argument of the maximization is the lower bound of the integral. The most promising approach that I've found seems to be "Quadrature on sparse grids":

http://sparse-grids.de/#Stata

which contains a Mata function nwspgr() in a zip file.  I've spent some time with this function and it's very accurate but I cannot figure out how to manipulate the integration bounds such that the lower bound is a variable that I can feed in from an optimizer such as Mata's optimize.

--
------------------------------------------------------------------------
Bob Hammond
Department of Economics
North Carolina State University
Office: (919) 513-2871
Fax: (919) 515-7873
http://www4.ncsu.edu/~rghammon/

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index