4310.sagews

def solve_for_weights_implicit():
    var("sigma,r,dt")
    mean=dt*(r-sigma^2/2);
    stddev=sqrt(dt)*sigma;
    var("lam,pUp,pDown");
    pMid=1-pUp-pDown;
    dx=dt*lam;
    assume(lam>0);
    m1= mean==pUp*dx-pDown*dx;
    m2= stddev^2== pUp*(dx-mean)^2+pMid*(-mean)^2+pDown*(-dx-mean)^2;
    sols=solve([m1,m2],[pUp,pDown])[0];
    u=sols[0].rhs().simplify_full();
    d=sols[1].rhs().simplify_full();
    m=(1-u-d).simplify_full()
    return u,m,d

class calc_weights:
    def __init__(self,_r,_sigma):
        u,m,d=solve_for_weights_implicit();
        discount=exp(-_r*dt);
        self.dx=(dt*lam).function(dt,lam);
        self.pu=u.subs(r=_r).subs(r=_r,sigma=_sigma).simplify_full().function(dt,lam);
        self.pm=m.subs(r=_r).subs(r=_r,sigma=_sigma).simplify_full().function(dt,lam);
        self.pd=d.subs(r=_r).subs(r=_r,sigma=_sigma).simplify_full().function(dt,lam);
        self.wu=(discount*u).subs(r=_r).subs(r=_r,sigma=_sigma).simplify_full().function(dt,lam);
        self.wm=(discount*m).subs(r=_r).subs(r=_r,sigma=_sigma).simplify_full().function(dt,lam);
        self.wd=(discount*d).subs(r=_r).subs(r=_r,sigma=_sigma).simplify_full().function(dt,lam);
 
    def __call__(self,_dt,_lam):
        return (self.wu(_dt,_lam),self.wm(_dt,_lam),self.wd(_dt,_lam));
        
    def mean(self,dt,lam):
        dx=self.dx(dt,lam);
        return (self.pu(dt,lam)*dx-self.pd(dt,lam)*dx)*dt;
        
    def stddev(self,dt,lam):
        dx=self.dx(dt,lam);
        mean=self.mean(dt,lam);
        return sqrt((self.pu(dt,lam)*(dx-mean)^2+self.pm(dt,lam)*(-mean)^2+self.pd(dt,lam)*(-dx-mean)^2));
    
    def skewness(self,dt,lam):
        dx=self.dx(dt,lam);
        mean=self.mean(dt,lam);
        stddev=self.stddev(dt,lam);
        tmp=(self.pu(dt,lam)*(dx-mean)^3+self.pm(dt,lam)*(-mean)^3+self.pd(dt,lam)*(-dx-mean)^3);
        return tmp / (stddev^(3/2));
    
    def kurtosis(self,dt,lam):
        dx=self.dx(dt,lam);
        mean=self.mean(dt,lam);
        stddev=self.stddev(dt,lam);
        tmp=(self.pu(dt,lam)*(dx-mean)^4+self.pm(dt,lam)*(-mean)^4+self.pd(dt,lam)*(-dx-mean)^4);
        return tmp / stddev^4;
        
    def metric1(self,dt,lam):
        eps_u=self.wu(dt,lam)-0.25;
        eps_m=self.wm(dt,lam)-0.5;
        eps_d=self.wd(dt,lam)-0.25;
        return eps_u^2+eps_d^2+eps_m^2;
        
    def metric2(self,dt,lam):
        eps_k=(self.kurtosis(dt,lam)-3)/3;
        eps_s=self.skewness(dt,lam);
        # Punishment factor to stop probs going negative
        err=max(0,-self.pu(dt,lam),-self.pm(dt,lam),-self.pd(dt,lam));
        return eps_k^2+eps_s^2+err^2*100000;

cw=calc_weights(0.02,0.1)

# Minimise for probabilities near (0.25,0.5,0.25)
vv=minimize(lambda x:cw.metric1(x[0],x[1]), [1,1]);
print(cw(vv[0],vv[1]))
print(cw.skewness(vv[0],vv[1]))
print(cw.kurtosis(vv[0],vv[1]))
tt=cw(vv[0],vv[1]);
print(log(abs(tt[0]-0.25))/log(2.0),log(abs(tt[1]-0.5))/log(2.0),log(abs(tt[2]-0.25))/log(2.0))

Warning: Maximum number of function evaluations has been exceeded. (0.258111879623, 0.500358951606, 0.241007716522) 9.56748169242e-05 2.00248387709 (-6.9457480411, -11.4439230258, -6.79709676795)

# Minimise for best skewness and kurtosis
vv=minimize(lambda x:cw.metric2(x[0],x[1]), [1,1]);
print(cw(vv[0],vv[1]))
print(cw.skewness(vv[0],vv[1]))
print(cw.kurtosis(vv[0],vv[1]))
tt=cw(vv[0],vv[1]);
print(log(abs(tt[0]-0.25))/log(2.0),log(abs(tt[1]-0.5))/log(2.0),log(abs(tt[2]-0.25))/log(2.0))

Optimization terminated successfully. Current function value: 0.000000 Iterations: 91 Function evaluations: 170 (0.120858032394, 0.653618430045, 0.205401045097) 1.48341044032e-07 3.00000128096 (-2.95297018164, -2.70257678417, -4.48684628623)