You want to expand it, cancelling quite often is done automatically

(diff(g(r),r)*f(phi)*r+r^2*diff(diff(g(r),r),r)*f(phi)+
  g(r)*diff(diff(f(phi),phi),phi)+h^2*g(r)*f(phi)*r^2)/g(r)/f(phi) = 0;
expand(%);

  1/g(r)*diff(g(r),r)*r + 
  1/g(r)*r^2*diff(diff(g(r),r),r) + 
  1/f(phi)*diff(diff(f(phi),phi),phi)+h^2*r^2 = 0

Though I usually do not enter such discussions: I agree, it is a nice example. It shows the deep understanding (of marketing? or more serious?) what 'academic' users need: nice pictures.

And I am not really sorry for that cynical remark.

At ~ 46 there is a pole, try simple plotting: (10088-10088*exp(-x))/(72*x+76)/((1-exp(-x))/x+10088/3/(72*x+76)-1);  plot(%, x=-1 .. 50);

Why you want that to be +1,0,-1 ? If those stand for stocks in an index,
say STOXX50, and stay for selecting, then usually one would do according
to their (relative) weight in the index.

The other is your (co)variance matrix. In which way ever you get it and
whatever it regards (dividends and tax as well?): this is not stable and
is always changing, so you recalculate and rebalance often.

One way to attack the NP complete problem would be to divide into smaller
parts by using sub-portfolios, say by clustering into sectors: stocks in
the same sector (for example Banks, Automotive, Health ...) usually are
highly correlated (mostly in a time stable way).

You will have ~ 7 clusters with ~ 7 stocks in it.

Then first you solve for one sector and then against your whole portfolio
(where you need fractions, that's why I asked for '+-1').

May be you have your own ideas, how to reduce the dimension(s).

This certainly is not a global optimum. But it has an economic rational
and it might come close to one - note, that no software will give you a
real global optimum.


However I am sure that there literature in Finance/Economics which does
cover that much more sophisticated.

Not sure, whether the package covers discrete problems as well.

Have you tried the brute way to set them as -1 , 0 , 1 before calling (ok, that would be 9 cases ..)?

Maple is not a reference and I do not want to step to such 'Standard' sheets,
sorry for my impatients. Somewhere s/he has to bring in the variable and for
that I would have to follow his coding (which is a mess in 2D and he also is
working with indexed variables). And autoexecution (which I hate for security
reasons). And 4 MB, while zipped it would be much smaller.

Here you go without obfuscation:

BSCall:= (s,e,t,r,v) ->
s*(1/2+1/2*erf(1/2*((ln(s/e)+r*t)/v/t^(1/2)+1/2*v*t^(1/2))*2^(1/2)))-
e*exp(-r*t)*(1/2+1/2*erf(1/2*((ln(s/e)+r*t)/v/t^(1/2)-1/2*v*t^(1/2))*2^(1/2)));

# log-normal Distribution
f:= (x,mu,sigma) -> exp(-((ln(x)-mu)/sigma)^2/2)/ (x*sigma*sqrt(2*Pi));

# un-discounted risk neutral density
RND:= (s,e,t,r,v) -> 1/2*exp(-1/2*((ln(s/e)+r*t)/v/t^(1/2)-
  1/2*v*t^(1/2))^2)*2^(1/2)/Pi^(1/2)/e/v/t^(1/2)

    s = spot price asset, current price of the underlying
    e = excercise price, strike
    t = maturity time, remaining time in years (annualized time)
    r = interest rate (continous annualized), 0.04 is 4 percent
    v = volatility (annualized), 0.20 is 20 percent
    p = option price, current price for an option

'f(strike,-v^2/2,v)/RND(1,strike,1,0,v)/exp(-0*1)' = 1;

                                 true


# natural bounds
assume(0<s):assume(0<e):assume(0<t):assume(0<v): assume(r::real):

# pricing with the RND, here for a Call
'exp(-r*t)*Int(max(S1-e,0) * RND(s,S1,t,r,v),S1=0..infinity)':
'exp(-r*t)*Int((S1-e)*RND(s,S1,t,r,v),S1=e..infinity)':
subs(s=1,t=1,r=0, %);
combine(%);
Change(%, S1=exp(kappa), kappa);
value(%);
%= 'BSCall(1,e,1,0,v)';
is(%);

                                 true

# finally show, that spot=1, time=1, rates=0 can be achieved:
'BSCall(s,e,t,r,v)'='s*BSCall(1,exp(-r*t)*e/s,1,0,v*sqrt(t))'; 
is(simplify(expand(%)));

                                 true


Whatever, the sheet does not give the price:

tstData:={r=0.050, v=0.5, s=100,t=1,K=100};
eval(BSCall(s,K,t,r,v), tstData); 
evalf(%);

                             21.79260422

while the sheet gives ~ 0.4 in (2.3) which is absurd.

The author should have 'felt' this. 

And a price is not given by 1 single erf. And a 'delta' is not
given by an exponential (but 1 single erf), since it is between
0 and 1, 'achiving' that for extreme strikes at the wings.

Any questions, why this is not a reference?

BTW: I even can not enter a new strike and hit the evaluation
button at (2.3) to get a new result. 

But yes, it looks nice, really. Certainly reviewed by marketing.

I do not understand stochastic processes at all, but your problem is based
on notational conventions (though you do not say what is what, r usually is
for rates ... and you do not gives a reference for your assertion).

Anyway: spot is a scaling (dito rates and time), strike is the variable.
You have to look into your derivation, one uses ln(strike) and  the other,
BS, uses ln(1/strike).

exp(-1/2*(ln(1/strike)/v-1/2*v)^2) = exp(-1/2*(ln(strike)+1/2*v^2)^2/v^2)
is what you want to see, Maple simplifies the quotient to 1 for 0 < strike.

For this you set spot=1, rates=0, time=1 (but if you work with the lognormal
then you know how to reduce to this case).

This is not readable. You may either use plain text of your input or upload a text file containing it if it is longer or even better upload a Maple sheet containing it (where the classical *.mws would be best).

If you use the following you will see, that the expression is monstrous:

theFinal:=diff(lmc,t):
length(theFinal)/80/60.0:
`pages`=floor(%);

                             pages = 3431

That means: the graphical interface would have to process ~ 3500 pages (for your 6 variables).
Which kills your application (and is useless I would guess).

The discounted risk neutral density is differentiating the BS twice w.r.t. to
the strike (which is called the Breeden-Litzenberger formula), RND:

'diff(BSCall(s,e,t,r,v),e$2)'='exp(-r*t)*diffN(dTwo(s,e,t,r,v))/(e*v*sqrt(t))'

dTwo:= (s, e, t, r, v) -> dOne(s,e,t,r,v)-v*t^(1/2),
dOne:= (s, e, t, r, v) -> (ln(s/e)+r*t)/v/t^(1/2)+1/2*v*t^(1/2)

diffN = D(N), N = x -> 1/2+1/2*erf(1/2*x*2^(1/2)) the cumulative normal fct.

If you set r=0 you have the undiscounted RND and get the pricing for almost
any option by price(strike) = exp(-r*t)*Int( payoff*RND(e), e in Reals), if
you tell Maple assume(0<s,0<e,0<v,0<t,r::real) and may be by helping through
narrowing the range (according to the payoff).

If the payoff is piecewise linear the answer should be a combinations of some
error functions, erf (since Maple does not use N itself)

I think your payoff is just 1 for 130 ... 140 and zero else, and r=0.

ahem ... your last side remarks would rule out your question as such :-)

anyway: it is the expectation value of the payoff (an indicator function for 130 ... 140 and paying 1 money unit) against the density (let me know, if you need to see it), which Maple would give if being feeded

If I remember correctly: under the assumptions of log-normality the 2nd terms in the Black-Scholes pricing formula can be seen as the probability of a Call resp Put of ending in the money - which should answer your question (setting rates = 0 and assuming you mean that as volatility). Or google for range options or (combined) digitals and their pricing formulae.

PS to get a feeling why that is quite unlikely: 12% is something like volatility for currencies (a bit higher), 100 may be seen as parity of USD and Euro. So the question is like: exchange rates increase by at least 30% (and at most 40%) within 1/2 a year. Usually that would need some very strong macro economical reason (which would increase volatility as well).

The command above can be understood as ploting a local chart (but normed to coordinates), since for a (smooth) surface you locally can always write it as f = 0.

I have never heard that one wants it through an atlas may some special graphical libraries can do, as you need to control the projections and glueings.

Are you aware of any existing software showing that?

PS: try to write plain text, it is a bit more easy to read that HTML ... (write in an text editor and cut&paste or choose formating 'normal')

It may be that questions are not answered which simply need some look up in the handbook or online help ...

If the surface is given through an explicit equation z = f(x,y) then you can use plot3d(f(x,y), x= - 3 .. 3, y= - 3 .. 3, axes=boxed)

If you have an implicit form you have to look up the help for the 3d version (do not have Maple here)