N:5 $ display(N)$

/* The Variance Gamma literature uses theta, nu, sigma, and omega = ln(1-theta nu -sigma^2 nu/2)/nu as parameters. */
/* Here, we use theta -> q  ,   nu -> v  ,   sigma -> c . */
/* We also use the abbreviation qs := theta + sigma^2/2 . */

load(distrib)$ powerdisp:true$
start_time():=block(time0:elapsed_real_time())$
stop_time():=block(time1:elapsed_real_time(), elapsed_time:float(round(1000.*(time1-time0))/1000), display(elapsed_time))$
define(Phi(x),cdf_normal(x,0,1))$
define(phi(x),ratsimp(diff(Phi(x),x)))$
b(x,s):=%e^(x/2)*Phi(x/s+s/2)-%e^(-x/2)*Phi(x/s-s/2)$
blacksigma:expand(fullratsimp(taylor(sum(sigma[k]*v^k,k,0,N),v,0,N)))$ display(blacksigma)$
define(psi_gamma(g,t,v),v^(-t/v)*g^(t/v-1)*%e^(-g/v)/gamma(t/v)); /* Unused. Only shown here for reference and decoration. */
define(phi_gamma(u,t,v),(1-%i*u*v)^(-t/v));
define(phi_gamma_expansion(u,t,v),taylor(phi_gamma(u,t,v),v,0,N));
singular_psi_gamma_expansion:taylor(fullratsimp(phi_gamma_expansion(dg/%i,t,v)/%e^(dg*t)),dg,0,2*N)$ display(singular_psi_gamma_expansion)$
/* In omega, we only keep terms to O(N) since all higher terms would be removed later anyway.
   It is assumed that doing this reduces memory usage and increases speed.                          */
omega:expand(fullratsimp((taylor(log(1-qs*v)/v,v,0,N))))$ display(omega)$

print("Computing Taylor expansion in nu of normalized Black call option price with sigma = sigma_0 + sigma_1*nu + ....")$
start_time()$
b_expansion_in_v:taylor(b(x,blacksigma*sqrt(t)),v,0,N)$
stop_time()$
print("Computing derivatives of normalized Black call option price conditioned on gamma variate value g.")$
print("Note that all these terms depend themselves on nu via omega whence we will need to expand them in nu later on.")$
start_time()$
define(b_conditioned_on_g[0](g),fullratsimp(  b(x+omega*t+qs*g,c*sqrt(g)) * %e^((omega*t+qs*g)/2)  ))$
for j:1 thru 2*N do define(b_conditioned_on_g[j](g),fullratsimp(diff(b_conditioned_on_g[j-1](g),g)))$
stop_time()$
print("Defining the normalized Variance Gamma call option price from the singular series.")$
start_time()$
define( vg(v), fullratsimp( sum( coeff(singular_psi_gamma_expansion,dg,j) * b_conditioned_on_g[j](t) , j , 0 , 2*N ) ) )$
stop_time()$
print("Expanding the Variance Gamma call option series once more in nu because all of the conditioned terms depend on nu themselves.")$
/* We are building the Taylor series manually, instead of using taylor(), to avoid some quirky messages otherwise issued by Maxima's taylor() function. */
define(vg_derivative[0](v),vg(v))$
for j:1 thru N do block( start_time(), define(vg_derivative[j](v),fullratsimp(diff(vg_derivative[j-1](v),v))), stop_time() )$
start_time()$
vg_expansion_in_v:sum(vg_derivative[j](0)*v^j/j!,j,0,N)$
stop_time()$

/* The result below, of course, means that sigma[0] = c. This was obviously expected, but we do it here as a consistency check. */

expand(fullratsimp( coeff(b_expansion_in_v,v,0) )); expand(fullratsimp(  coeff(vg_expansion_in_v,v,0) ));

/* Et voila ! */
print("")$
print("Recall that:")$
display(blacksigma)$
print("      c   =   sigma")$
print("      v   =   nu")$
print("      qs  =   theta + sigma^2/2")$
solution[0]:sigma[0]=c$
for j:1 thru N do block(
 tmp[0]: fullratsimp(coeff(b_expansion_in_v,v,j)),
 for k:1 thru j do tmp[k]:fullratsimp(subst(rhs(solution[k-1]),sigma[k-1],tmp[k-1])),
 solution[j] : expand(solve( fullratsimp(coeff(vg_expansion_in_v,v,j)) = tmp[j] , sigma[j] )[1])
)$
for j:0 thru N do block(print(""),display(solution[j]))$

define(sigma_of_x(x),sum(rhs(solution[j])*v^j,j,0,N))$

sigma_atm:block([tmp:taylor(sigma_of_x(0),v,0,N)],sum(expand(coeff(tmp,v,j))*v^j,j,0,N))$ display(sigma_atm)$

dsigmadx_atm:block([tmp:taylor(subst(0,x,diff(sigma_of_x(x),x)),v,0,N)],sum(expand(coeff(tmp,v,j))*v^j,j,0,N))$ display(dsigmadx_atm)$

d2sigmadx2_atm:block([tmp:taylor(subst(0,x,diff(sigma_of_x(x),x,2)),v,0,N)],sum(expand(coeff(tmp,v,j))*v^j,j,0,N))$ display(d2sigmadx2_atm)$

dsigmadk_atm_over_k:block([tmp:taylor( - dsigmadx_atm ,v,0,N)],sum(expand(coeff(tmp,v,j))*v^j,j,0,N))$ display(dsigmadk_atm_over_k)$

d2sigmadk2_atm_over_k2:block([tmp:taylor( dsigmadx_atm + d2sigmadx2_atm ,v,0,N)],sum(expand(coeff(tmp,v,j))*v^j,j,0,N))$ display(d2sigmadk2_atm_over_k2)$