Maxima crashes on relatively simple integral

Question

I'm trying to Maxima-fy my Mathematica box options formula (https://github.com/barrycarter/bcapps/blob/master/box-option-value.m) but Maxima crashes on a fairly simple integration:

load(distrib); 
pdflp(x, p0, v, p1, p2, t1, t2) := pdf_normal(x,log(p0),sqrt(t1)*v); 
cdfmaxlp(x, p0, v, p1, p2, t1, t2) := 1-erf(x/(v*sqrt(t2-t1)/sqrt(2))); 

upandin(p0, v, p1, p2, t1, t2) :=  
 integrate( 
 float( 
 pdflp(x, p0, v, p1, p2, t1, t2)* 
 cdfmaxlp(log(p1)-x, p0, v, p1, p2, t1, t2) 
 ), 
 x, minf, log(p1)); 

Evaluating upandin w/ certain values crashes:

upandin(1, .15, 1.01, 1.02, 1/365.2425, 2/365.2425); 

rat: replaced -.00995033085316809 by -603/60601 = -.00995033085262619 

rat: replaced 2.718281828459045 by 23225/8544 = 2.718281835205993 

rat: replaced 8116.5 by 16233/2 = 8116.5 

rat: replaced 2.718281828459045 by 23225/8544 = 2.718281835205993 

rat: replaced -8116.5 by -16233/2 = -8116.5 

rat: replaced 1.0 by 1/1 = 1.0 

rat: replaced 1.792882852833688 by 4484/2501 = 1.792882846861255 

rat: replaced 180.1832400641081 by 126849/704 = 180.1832386363636 

rat: replaced 2.718281828459045 by 23225/8544 = 2.718281835205993 

rat: replaced -8116.5 by -16233/2 = -8116.5 

rat: replaced -1.0 by -1/1 = -1.0 

rat: replaced 1.792882852833688 by 4484/2501 = 1.792882846861255 

rat: replaced 180.1832400641081 by 126849/704 = 180.1832386363636 

rat: replaced 2.718281828459045 by 23225/8544 = 2.718281835205993 

rat: replaced -8116.5 by -16233/2 = -8116.5 

rat: replaced 1.0 by 1/1 = 1.0 

rat: replaced -1.0 by -1/1 = -1.0 
Maxima encountered a Lisp error: 

 The value 16090668801 is not of type FIXNUM. 

Without the float() in upandin, Maxima just leaves the integral in original form.

Can someone help? I thought converting Mathematica to Maxima would be easy, but now I'm not as sure.

The Mathematica version works fine:

pdflp[x_, p0_, v_, p1_, p2_, t1_, t2_] :=  
 PDF[NormalDistribution[Log[p0],Sqrt[t1]*v]][x] 

cdfmaxlp[x_, p0_, v_, p1_, p2_, t1_, t2_] := 1-Erf[x/(v*Sqrt[t2-t1]/Sqrt[2])]; 

(* NIntegrate below "equivalent" to Maximas float(); no closed form *) 

upandin[p0_, v_, p1_, p2_, t1_, t2_] :=  
 NIntegrate[pdflp[x, p0, v, p1, p2, t1, t2]* 
           cdfmaxlp[Log[p1]-x, p0, v, p1, p2, t1, t2], 
{x, -Infinity, Log[p1]}] 

upandin[1, .15, 1.01, 1.02, 1/365.2425, 2/365.2425] 

0.0998337 

EDIT: Is there any open source Mathematica-like program that WILL numerically approximate this function? I'd really like to release open source code to an open source platform.

Answer 1

Use quad_qagi to numerically approximate an integral over an infinite interval. ?? quad_ shows info about Quadpack functions.

load (distrib);
pdflp (x, p0, v, p1, p2, t1, t2) := pdf_normal (x, log(p0), sqrt(t1)*v); 
cdfmaxlp (x, p0, v, p1, p2, t1, t2) := 1 - erf(x/(v * sqrt(t2 - t1)/sqrt(2))); 

upandin (p0, v, p1, p2, t1, t2) := block ([integrand],
   integrand : pdflp (x, p0, v, p1, p2, t1, t2) * cdfmaxlp (log(p1) - x, p0, v, p1, p2, t1, t2),
   quad_qagi (integrand, x, minf, log(p1))); 

upandin (1, .15, 1.01, 1.02, 1/365.2425, 2/365.2425);
 => [.09983372557898755, 2.839204848435967E-10, 225, 0]

Sorry for the late reply. Leaving this here in case someone finds it by searching.

Answer 2

Maxima's 'integrate' function does symbolic, not numeric, integration. When it returns a noun form from an integral, that means it can't perform the (symbolic) integration. Changing the parameters of the expression from exact to floating (using 'float') won't change that.

I think what you're looking for is a numeric integration routine -- Maxima offers a variety of those, from the very basic romberg to a variety of Quadpack methods (try ?? quad for documentation).

    -s

PS As for "this crappy 'open source' stuff" -- what brought that on? You might want to look at the history of Macsyma/Maxima in the Wikipedia article for some perspective.

Answer 3

I know Maxima tries very hard to avoid floats, and I think that's what it's trying to do here, but I'm not enough of a Maxima guru to explain how to prevent it. Pretty much anything numerical can handle this, although you might have to break the interval or transform the integrand manually. Note that you say it's fairly simple, but it's awfully steep: for these parameters, the integrand is ~6*10^(-34) at 0.1 and ~3*10^(-206) at -0.1. That's enough of a range to give lots of naive integration algorithms fits.

Anyway, you can do it easily enough in Sage using tools from scipy and gsl behind the scenes:

import scipy.stats

def pdflp(x,p0,v,t1):
    return scipy.stats.norm(log(p0), sqrt(t1)*v).pdf(x)

def cdfmaxlp(x,v,t1,t2):
    return (1-erf(x/(v*sqrt(t2-t1)/sqrt(2.))))

def upandin(p0, v, p1, p2, t1, t2):
    integrand = lambda x: (pdflp(x,p0,v,t1) * cdfmaxlp(log(p1)-x,v,t1,t2))
    return numerical_integral(integrand, -Infinity, log(p1))

sage: upandin(1, .15, 1.01, 1.02, 1/365.2425, 2/365.2425)
(0.099833725578983457, 7.5174412058308382e-07)

or use mpmath's quad if you need arbitrary precision. [I had a guess at the "correct" value here, but since we don't have that much precision to start with, it's kind of silly.]

Answer 4

(I probably have no business answering this, but...)

Just a guess, but it seems that integrate wants to make the input exact again, and maybe is doing some difficult bignum computations involving rational arithmetic. It rationalize your approximate e (Euler number) so that means it could behave differently from integrate(0 with exact input.

Might want to check

http://eagle.cs.kent.edu/MAXIMA/maxima_21.html

or

http://www.delorie.com/gnu/docs/maxima/maxima_62.html

for dedicated numerical code e.g from Quadpack.

(Still wondering why I'm even trying to answer this. There must be Maxima expertise somewhere on Stack Overflow.)

Daniel Lichtblau Wolfram Research