# Gram-Charlier Processes and Applications to Option Pricing

Authors: Jean-Pierre Chateau; Daniel Dufresne
Journal: Journal of Probability and Statistics
 (2017)

---

## Abstract

A Gram-Charlier distribution has a density that is a polynomial times a normal density. For option pricing this retains the tractability of the normal distribution while allowing nonzero skewness and excess kurtosis. Properties of the Gram-Charlier distributions are derived, leading to the definition of a process with independent Gram-Charlier increments, as well as formulas for option prices and their sensitivities. A procedure for simulating Gram-Charlier distributions and processes is given. Numerical illustrations show the effect of skewness and kurtosis on option prices.


---

## Body

## 1. Introduction

Gram-Charlier series are expansions of the form(1)fx=ϕa,b;x1+c1He1x-ab+c2He2x-ab+⋯,where (2)ϕa,b;x=1b2πe-1/2b2x-a2,x∈R,is the usual normal density and Hek is the Hermite polynomial of order k. The expression within square brackets in (1) is an orthogonal polynomial expansion for the ratio f(x)/ϕ(x); given an arbitrary function the expansion may or may not converge to the true value of f(x). In this paper we focus on the properties of the Gram-Charlier distributions, obtained by truncating the series after a finite number of terms. What is obtained is a family of distributions parametrized by a,b,c0,…,cN, as is explained in detail below.This paper has three main goals: (1) define and study the properties of the family of Gram-Charlier distributions; (2) define a Gram-Charlier process and derive its basic properties; (3) apply those to European options. The formulas we give for European option prices and Greeks apply to Gram-Charlier distributions of any order, and we use four- and six-parameter Gram-Charlier distributions in our examples. Two numerical illustrations show how option prices are affected by the skewness and kurtosis of returns. This paper can a reference for those using Gram-Charlier distributions in option pricing but also in statistics.Most previous applications to option pricing have assumed thatc1=c2=0. We believe this restriction is not necessary in a general theory of Gram-Charlier distributions, as there may well be situations where the extra degrees of freedom given by c1 and c2 will be useful. We discuss it in detail at the end of Section 2.It has been observed that option prices have nonconstant implied volatilities, meaning that log returns do not have a normal distribution under the risk-neutral measure. There is a wide literature on modelling log returns to fit observed option prices, the main alternatives to Brownian motion being stochastic volatility models (where the parameterσ in Black-Scholes is replaced with a continuous-time stochastic process), GARCH time series, and Lévy processes. Gram-Charlier distributions are mathematically simpler than the models just mentioned, while allowing a better fit to data than the normal distribution. Several authors have used Gram-Charlier distributions in option pricing, as a model for log returns, among others [1–10]. The majority of previous authors assumed a density of the form (3)ϕx1+s¯6He3x+k¯24He4x,ϕx=ϕ0,1,x=12πe-x2/2,for the normalized log return. (In our notation this is a GC(0,1;0,0,s¯/6,k¯/24) distribution; see Section 2.) The notation emphasizes that in this case the coefficient of He3(x) turns out to be Pearson’s skewness coefficient s¯ divided by 6, and the coefficient of He4(x) the excess kurtosis coefficient k¯ divided by 24. The distribution of log returns is then a four-parameter Gram-Charlier distribution (since there are two other parameters, the mean and variance, besides s¯ and k¯). This distribution allows nonzero skewness and excess kurtosis, unlike the normal distribution found in Black-Scholes. In (3) the parameters (s¯,k¯) are restricted to a specific region (see Figure 3), because outside that region the function in (3) becomes negative for some values of x (see Section 2.1).Jurczenko et al. [8] specify the martingale restriction (4)S0=e-rTEQSTthat the four-parameter Gram-Charlier density must satisfy in pricing options (previous authors had not taken it into account). The martingale condition for general Gram-Charlier distributions is described in Section 5. Our Gram-Charlier distribution with parameters a,b,c1,…,cN, denoted by GC(a,b;c1,…,cN), has density (5)ϕa,b,x1+c1He1x-ab+⋯+cNHeNx-ab.As stated above we do not set c1=c2=0, and N can be any even positive integer. The question whether (5) is nonnegative for all x is an important one. Several authors have disregarded this issue and, in fact, some have come up with parameters that do not yield a true probability density function. One might argue that this is the price to pay for truncating an infinite Gram-Charlier series. However, if the same log return distribution is used to price many options, then a true probability density function is the only safe choice, because otherwise there might be inconsistencies among option prices. For instance, if the function used as density is negative over the interval (α,β), then a digital option that pays off only when the log return is in that interval will have a negative price. In this paper we consider Gram-Charlierdistributions, notexpansions, and insist that the densities integrate to one and be nonnegative. Our goal is to define a family of proper probability distributions; nevertheless, our formulas do apply to truncated Gram-Charlier expansions as well.The paper by León et al. [11] presents an alternative to the general Gram-Charlier distributions we study in this paper. Those authors consider the subclass of Gram-Charlier distributions consisting of densities(6)fXx=ϕa,b,xpx,where the polynomial p(x) is the square of another polynomial, p(x)=q(x)2. This has the obvious advantage that the nonnegativity restriction on fX(·) is automatically satisfied. We discuss that subclass of “squared” Gram-Charlier distributions in the Conclusion.Almost all previous authors have used Gram-Charlier distributed log returns over a single time period. This has an obvious downside, in that it becomes tricky, if not impossible, to preserve consistency between the prices of options with different maturities. Section3 shows that a Lévy process with Gram-Charlier increments does not exist; however, it also shows that the sum of independent Gram-Charlier distributed variables also has a Gram-Charlier distribution. This opens the way for multiperiod Gram-Charlier option pricing, using a discrete-time random walk model for which the log return over any period has a Gram-Charlier distribution. The Gram-Charlier distribution of the multiperiod return has a larger number of parameters, though the model is still simpler than almost any (if not all) alternative stochastic volatility models. There is no problem computationally, since we give explicit formulas for options under Gram-Charlier distributions with an arbitrary number of parameters.The layout of the paper is as follows. In Section2, we extend the study of Gram-Charlier distributions to all possible polynomials p(·) and derive their properties (moments, cumulants, moment determinacy, properties of the set of valid parameters, tail, and so on). Some of the formulas and properties in Theorem 2 appear to be new. In Section 3 we show that there is no Lévy process with Gram-Charlier distributed increments, apart from Brownian motion, and we define a discrete-time process with independent Gram-Charlier increments that is suitable for option pricing. In Section 4 we show that the log Gram-Charlier distribution is not determined by its moments, just like the lognormal. Next, Section 5 gives formulas for European call and put prices when the log-price returns of the underlying have a general Gram-Charlier distribution; the previous literature mostly considered the GC(a,b;0,0,c3,c4) family and the squared Gram-Charlier distributions (an exception is [10], where it is assumed that c1=c2=0). In particular, we derive a change of measure formula that extends the Cameron-Martin formula for the normal distribution; the latter is used in pricing European options in the Black-Scholes model. We also derive formulas for the sensitivities (Greeks) of those option prices with respect to all parameters. A technique for simulating Gram-Charlier distributions is described. Parts (c), (d), (i), (j) (k), and (l) of Theorem 2 and part (b) of Theorem 3 and Theorems 6, 7, 9, and 11 appear to be new. Theorem 3(a) is for the first time formulated for general Gram-Charlier distributions. Theorems 3(a), 4, and 8 are given, more or less explicitly, for the subclass of squared Gram-Charlier distributions in [11], and some of the Greeks in Theorem 9 had been calculated for the GC(a,b;0,0,c3,c4) distribution by previous authors.In Section6 we give two applications that show how option prices depend on skewness and kurtosis of the log returns (this of course cannot be done in the classical Black-Scholes setting, while it appears quite complicated to do so in stochastic volatility or Lévy driven models). The first example is equity indexed annuities (EIAs in the sequel) premium options. The pricing of EIAs has been studied by many authors, including Hardy [12], Gaillardetz and Lin [13], and Boyle and Tian [14]. The second example is lookback options, which illustrates the use of the simulation method in Section 2.5. Lookback options have been studied by numerous authors, but there is no closed form formula for the price of discretely monitored lookback options within the Black-Scholes model. We refer the reader to Kou [15]. In those examples all parameters are estimated by maximum likelihood, with the range of the parameters restricted to where they correspond to a true probability distribution. Moment estimation is of course possible when c1=c2=0, but it is not trustworthy, because of the range restriction on the parameters c3,c4,…; for instance, the empirical skewness and excess kurtosis have a positive probability of falling outside the feasible region (see Section 6). Theorem 2(f) shows that the first k moments of the Gram-Charlier distribution depend on more than k parameters, which probably rules out moment estimation, unless one a priori fixes two of the parameters (which again is not trustworthy). We apply maximum likelihood to the six-parameter GC(a,b;c1,c2,c3,c4) in Section 6, but we do realize that maximum likelihood estimation for higher order Gram-Charlier distributions poses computational problems, which are an interesting avenue for further research. (A referee pointed out that fixing a and b would make estimation easier; our guess is that if one wishes to simplify maximum likelihood estimation then fixing a and b might be a better idea than fixing c1 and c2, though we have not looked at this in any depth. A two-step process might be imagined, whereby the data first give information about a and b, and then maximum likelihood is applied using that preliminary information. This would help by constraining the optimization to a smaller region, while agreeing with the intuitive idea that a and b are location and scale parameters, resp.)Notation 1.
The normal density function is denoted as(7)ϕa,b,x=1b2πe-1/2b2x-a2,ϕ0,1,x=ϕx=e-x2/22π,and its distribution function is (8)Φx=∫-∞xϕydy.Two equivalent versions of the Hermite polynomials may be found in the literature: for n=0,1,2…, (9)Hnx=-1nex2dndxne-x2Henx=-1nex2/2dndxne-x2/2.The first one is common in mathematics and physics, but in probability and statistics there is an obvious advantage in using the second one. (The conversion formula is Hen(x)=2-n/2Hn(x/2).) The first few Hermite polynomials are (10)He0x=1,He1x=x,He2x=x2-1,He3x=x3-3x,He4x=x4-6x2+3,He5x=x5-10x3+15x,He6x=x6-15x4+45x2-15.


## 2. Gram-Charlier Distributions

For a fixedN, consider the class of distributions that have a pdf of the form(11)fx=ϕx∑k=0NckHekx,x∈R,with cN≠0. Noting that the leading term of Hek(x) is xk, we conclude that N must necessarily be even, because if N were odd then the polynomial that multiplies ϕ(x) would take negative values for some x. For the same reason cN cannot be negative.Definition 1.
Leta∈R, b>0, ck∈R, c0=1, and N∈{0,2,4,…}. We write Y~GC(a,b;c1,…,cN) (or Y~GC(a,b;c~)) if the variable (Y-a)/b has probability density function(12)ϕx∑k=0NckHekx.This will be called aGram-Charlier distribution with parameters a,b,c~, with c~=(c1,…,cN). The largest N such that cN>0 is called theorder of the Gram-Charlier distribution. The normal distribution with mean a and variance b2 is a GC(a,b;(0,…0)) (or GC(a,b;-)) with order 0.The class of Gram-Charlier distributions just defined includes all distributions with density(13)1bϕy-abpy,where p(y) is a polynomial of degree N, since p(y) can be rewritten as a combination of (14)Heky-ab,k=0,1,…,N.The condition c0=1 ensures that function (12) integrates to one, since (15)∫-∞∞ϕxHekxdx=-1k∫-∞∞dkdxkϕxdx=0,k=1,2,….There are no simple conditions that ensure that a polynomial remains nonnegative everywhere, though in some cases precise conditions on c~ are known; see below. If a vector c~ leads to a true Gram-Charlier pdf, then we will say that c~ isvalid.Generating functions are convenient when dealing with orthogonal polynomials. One is(16)w1t,x=∑n=0∞tnn!Henx=etx-t2/2.Another one is (17)w2t,u,x,y=∑k=0∞∑n=0∞tkk!unn!∫-∞∞e-x2/2HekxHenx+ydx=2πeut+y.Letting y=0 leads to (18)12π∫-∞∞e-x2/2HekxHenxdx=∂ktk∂nunetut=u=0=0ifk≠nn!ifk=n.This proves the orthogonality of the Hermite polynomials and gives us the value of (19)12π∫-∞∞e-x2/2Hen2xdx,which is essential in deriving Gram-Charlier series. Since He0=1 this also implies(20)∫-∞∞ϕxHekxdx=0,k=1,2,….Another formula is the Laplace transform of ϕ(x)Hen(x), which may be found by integrating by parts n times:(21)∫-∞∞etxϕxHekxdx=tket2/2,k=0,1,….Let us calculate the moments of a distribution with density (11). First consider(22)∫-∞∞xnϕxHekxdx=-1k∫-∞∞xndkdxkϕxdx.Integrating by parts repeatedly yields(23)∫-∞∞xnϕxHekxdx=n!2n-k/2n-k/2!ifn-kis  even  and  nonnegative0ifn-kis  odd  and  nonnegative.Hence, the nth moment of the distribution in (11) is (24)∑k=0N∧nckn!2n-k/2n-k/2!1n-keven.This says in particular that the parameter ck only affects the moments of order k and higher of the GC(a,b;c~) distribution (see part (i) of Theorem 2). This is confirmed by the moment generating function, which may be found from (21): (25)∫-∞∞etxfxdx=et2/2∑k=0Ncktk.It can be checked that differentiating this expression n times and setting t=0 give the same expression found for the nth moment of (11).Theorem 2.
SupposeY~GC(a,b;c~), c~∈RN with b>0, c0=1, and cN>0. The order N of the distribution is necessarily even.
(a)(26)EY-an=bn∑k=0N∧nckn!2n-k/2n-k/2!1n-keven.
(b)(27)EetY=eat+b2t2/2∑k=0Nbkcktk,t∈R.
(c) The representation of theGC(a,b;c~) distribution in terms of the parameters a, b, and c~ is unique.
(d) All Gram-Charlier distributions are determined by their moments.
(e) The set of validc~ in RN includes the origin, is not reduced to a single point, and is convex.
(f) The first six moments of theGC(a,b;c~) distribution are (28)m1=a+bc1,m2=a2+b2+2abc1+b2c2,m3=a3+3ab2+3bc1a2+b2+6ab2c2+b3c3,m4=a4+4a3bc1+12a2b2c2+6a2b2+12ab3c1+24ab3c3+12b4c2+24b4c4+3b4,m5=a5+5a4bc1+20a3b2c2+10a3b2+30a2b3c1+60a2b3c3+60ab4c2+120ab4c4+15ab4+15b5c1+60b5c3+120b5c5,m6=a6+6a5bc1+30a4b2c2+15a4b2+60a3b3c1+120a3b3c3+180a2b4c2+360a2b4c4+45a2b4+90ab5c1+360ab5c3+720ab5c5+90b6c2+360b6c4+720b6c6+15b6.
(g) The first six cumulants of theGC(a,b;c~) distribution are (29)κ1=a+bc1,κ2=b21-c12+2c2,κ3=2b3c13-3c1c2+3c3,κ4=-6b4c14-4c12c2+2c22+4c1c3-4c4,κ5=24b5c15-5c13c2+5c1c22+5c12c3-5c2c3-5c1c4+5c5,κ6=-120b6c16-6c14c2+9c12c22-2c23+6c13c3-12c1c2c3+3c32-6c12c4+6c2c4+6c1c5-6c6.
(h) The following hold for anyGC(a,b;c~) distribution:(30)mean:a+bc1variance:b21-c12+2c2skewnesscoefficient:2c13-3c1c2+3c31-c12+2c23/2excesskurtosiscoefficient:-6c14-4c12c2+2c22+4c1c3-4c41-c12+2c22.
(i) SupposeX~GC(a,b;c~X), Y~GC(a,b;c~Y). Then the first K moments of X and Y are the same; that is,(31)EXj=EYj,j=1,…,K,if and only if cjX=cjY, j=1,…,K. This implies that if c1 up to c4 are equal to 0 then the distribution has zero skewness and excess kurtosis; hence, this shows how to construct an infinite number of distributions that share this property with the normal.
(j) SupposeX~GC(a,b;c~X). Then (32)a=EX⟺c1X=0,b2=EX-a2⟺c2X=0,a=EX,b2=VarX⟺c1X=c2X=0.When c1X=c2X=0 the skewness and excess kurtosis coefficients of X are 6c3X and 24c4X, respectively, for any N=0,2,4,….
(k) IfX~GC(a,b;c1,…,cN) and q is a constant then Y=qX~GC(aq,b|q|;c1′,…,cN′), where ck′=(sign⁡(q))kck, k≥1. In particular, -X~GC(-a,b;-c1,…,(-1)N-1cN-1,cN).
(l) The law of the square ofX~GC(0,1;c1,…,cN) is a combination of chi-square distributions with 1,3,…,N+1 degrees of freedom that has density (33)gy=e-y/22πy∑j=0n/2αjyj,where∑j=0n/2αjyj=∑j=0n/2c2jHe2jy.Proof.
Parts (a) and (b) were proved above, and (c) follows directly from (b). To prove (d) it is sufficient to note the existence ofEetY for t in an open neighbourhood of s=0. For (e), if c~=0~∈RN then the distribution is the standard normal, and this is a Gram-Charlier distribution. Since cN>0 and N is even, the polynomial(34)∑k=1NHekxtends to ∞ when |x| tends to ∞. Hence, there is ϵ0>0 such that(35)1+ϵ∑k=0NHekx>0,x∈R,0≤ϵ≤ϵ0(recall that c0He0(x)=1 for all Gram-Charlier distributions). The set of valid vectors c~ thus includes (ϵ,…,ϵ)∈RN for each 0≤ϵ≤ϵ0. If c~(j),j=1,2, are valid then pc~(1)+(1-p)c~(2) is also valid, for any p∈(0,1).
Part (f) is found by expanding the mgf in (b) as a series ins around the origin. Part (g) is found by expanding (36)t⟼logeat+b2t2/2∑k=0Nbkcktk.The formulas in (h) follow (f) and (g).
For part (i), it is sufficient to consider the casea=0, b=1 only. Write(37)ψt=et2/2,qXt=∑k=0NckXtk,qYt=∑k=0NckYtkand calculate the moments of X and Y by successively differentiating the mgf’s of X and Y (note that c0X=qX(0)=c0Y=qY(0)=1). The first K moments of X and Y are the same if and only if(38)ψ′0qX0+ψ0qX′0=ψ′0qY0+ψ0qY′0,⋮∑j=0KKjψj0qXK-j0=∑j=0KKjψj0qYK-j0.Suppose that cjX=cjY, j=1,…,K. Then qX(j)(0)=qY(j)(0), j=1,…,K, and thus the first K moments of X and Y are the same. Conversely, suppose EXj=EYj,j=1,…,K. Then the first identity above implies that c1X=qX′(0)=qY′(0)=c1Y, since ψ(0) is not zero. The second identity implies that c2X=c2Y, and so on, up to cKX=cKY.
Turning to (j), the first equivalence(39)a=EX⟺c1X=0is an immediate consequence of property (i) with K=1. To prove the second equivalence, suppose X~GC(0,1;c~), and let Y~GC(0,1;-)=N(0,1). Then (40)EX2=ψ′′0qX0+2ψ′0qX′0+ψ0qX′′0,EY2=ψ′′0qY0+2ψ′0qY′0+ψ0qY′′0.Here qX(0)=qY(0)=1 and ψ′(0)=0; hence, EX2=EY2 if and only if qX′′(0)=qY′′(0). The last equality is c2X=c2Y. For X~GC(a,b;c~X), with a and b being arbitrary, this means that (41)Ex-ab2=1⟺c2X=0.For (k), observe that if q≠0 then the mgf of Y is (42)eaqs+1/2b2q2s2∑k=0Nckbkqksk=eaqs+1/2b2q2s2∑k=0Nckqqkbqksk.Finally, turn to (l). Routine calculations show that the density of the square is (43)12yfXy+fX-y=ϕy2y∑k=0NckHeky+He-y.The Hermite polynomials of odd order are odd functions and so disappear from that expression, while the even order Hermite polynomials are even functions.WhenN=0 the GC(a,b;c~) distribution is the normal distribution with mean a and variance b2. However, part (b) of the theorem says that when N>0 the parameters a and b2 are not necessarily the mean and variance of the distribution. Simple calculations lead to the following result.Theorem 3.
(a) IfX~GC(0,1;c1,…,cN) then (44)PX≤x=Φx-ϕx∑k=1NckHek-1x,PX>x=Φ-x+ϕx∑k=1NckHek-1x.
(b) The tails of theGC(a,b;c1,…,cN) distribution are (45)PX>x~cNx-abN-1ϕx-abasx⟶∞,PX<x~cNx-abN-1ϕx-abasx⟶-∞.The tails of the Gram-Charlier distributions are thicker than those of the normal distribution but are still “thin” because they are in the limit smaller than any exponential functionexp⁡(-αx).


### 2.1. TheGC(a,b;0,0,c3,c4) Family

Here the exact region for the(c3,c4) that lead to a true probability distribution has been found. This goes back to Barton and Dennis [16], but a more detailed explanation is given in Jondeau and Rockinger [6]. The region is shown in Figure 3 (use the correspondence s¯=6c3, k¯=24c4 from Theorem 2(h)). An important fact about this region is that it is not rectangular; the possible excess kurtosis values depend on skewness, and conversely.


### 2.2.GC(a,b;c~) Distributions of Order 4 and Higher

TheGC(a,b;c1,c2,c3,c4) family has six parameters, rather than four, and thus has more degrees of freedom in fitting data; to the authors’ knowledge the general six-parameter Gram-Charlier distribution has been used in financial applications by León et al. [11] only (those authors use the subfamily consisting of polynomials p(·) that are squares of some second-degree polynomial). Schlögl [10] fits the six-and eight-parameter families GC(a,b;0,0,c3,c4,c5,c6) and GC(a,b;0,0,c3,c4,c5,c6,c7,c8) to data.The set of(c1,c2,c3,c4) that yield true probability distributions has not been identified, but it is possible to fit the six parameters a,b,c1,c2,c3,c4 by maximum likelihood.


### 2.3. Why We Do Not Assume Thatc1=c2=0

Almost all previous authors have assumed thatc1=c2=0, because they used normalized data:(46)Y=X-EXVarX(see part (j) of Theorem 2) as explained below. In this section we explain why it is important not to restrict Gram-Charlier distributions or series in that way. The first reason for not doing so is that enlarging the parameter space can only be a good thing. The second one is that in fitting those distributions to data there may be very real advantages in letting c1 and c2 be different from zero. The only downside of letting c1 and c2 take nonzero values, and it is of no real importance, is that c3 and c4 lose their simple relationship with skewness and excess kurtosis (see part (h) of Theorem 2). A third reason is that after an exponential change of measure a Gram-Charlier distribution will rarely have c1=c2=0; see Section 2.4.We now show, using an example that can be worked out explicitly, that it is not always best to use normalized data when fitting Gram-Charlier distributions, because choosing another affine transformation of the data may well yield a much better fit.The “standard” Gram-Charlier expansion for a functionf(·) is(47)ϕx∑k=0∞ckHek,where (48)ck=1k!∫-∞∞Hekxfxdx=1k!EHekX.A classical result about Hermite series, proved by Cramér [17], is that sufficient conditions for the Gram-Charlier expansion (47) to converge to (49)12fx-+fx+for all x∈R are that (i) f(·) has finite variation in every bounded interval, and (ii) satisfies (50)∫-∞∞ex2/4fxdx<∞.If X has density f(·), then the last condition is EeX2/4<∞. This condition cannot be improved upon, in the sense that there are cases where the Gram-Charlier series defined above diverges, although EeρX2<∞ for all ρ<1/4 (it will be shown below that one such case is the normal distribution with mean 0 and variance 2).Let us first calculate(51)hq,k=EHekX,X~N0,q2.Using the generating function (16) we find (52)hq,k=0ifkis  oddk!2k/2k/2!q2-1k/2ifkis  even.The Gram-Charlier expansion corresponding to the N(0,q2) density is thus(53)ϕx∑k=0∞hq,kk!Hekx=ϕx∑n=0∞12nn!q2-1nHe2nx.It is possible to determine whether this converges or not when x=0; the Gram-Charlier series for f(0) is (54)12π∑n=0∞-1n2n!22nn!2q2-1n.From Stirling’s formula, (55)2n!22nn!2~1πnas n tends to infinity, which implies that there are three possibilities: expansion (54) (i) converges absolutely if 0<q2<2, (ii) converges simply if q2=2, and diverges if q2>2. (The preceding calculations are from Cramér [17].)Let us now consider a random variableXp, 0≤p≤1, with density (56)gpx=pϕx+1-p1pϕxp.In words, Xp has a N(0,1) distribution with probability p and a N(0,p2) distribution with probability 1-p. Define the normalized variable (57)Yp=XpVarXp=Xpσp,σp2=VarXp=p+p2-p3.The coefficients cp,k of the Gram-Charlier expansion for Yp are (58)cp,k=1k!EHekYp=1k!ph1/σp,k+1-php/σp,k.We have shown that (53) diverges if q2>2. Here we see that cp,k involves hq,k with q=1/σp and q=p/σp; in the second case there is no problem, as (p/σp)2<2 for all p∈(0,1]. In the first case, (59)1σp2>2⟺p+p2-p3<12.This condition is satisfied for p smaller than p∗≈0.403, and so the Gram-Charlier expansion for Yp diverges for all p<p∗. This is a sad state of affair: a series designed to work for distributions that are “close to the normal” that fails for combinations of two normal densities!There is an easy solution: use a different scaling forXp. Rather than multiplying Xp by 1/σp, use a factor r such that the series converges. In this example we know that the series converges if and only if r<2. Say we choose p=.3; this implies (60)1VarXp=10.3+0.32-0.33=2.75482>2.Figure 1 shows how the true density for X0.3 is approximated by order 10 and 20 Gram-Charlier truncated series if r=1/σ0.3≐1.66 is chosen; the graph is not surprising, the infinite Gram-Charlier expansion diverges. Figure 2 shows the same except that r is set to 1.4. The latter is smaller than 2, so the infinite expansion converges. In this example the density to be approximated is symmetrical about 0, so c1 is always 0, while (61)q=1σ0.3⟹c2=0;q=1.4⟹c2=-0.14426.Figure 1
Gram-Charlier approximations to the densityg0.3(x) (solid line) when q=1/σ0.3≐1.66, resulting in c1=c2=0. The curve with long dashes is the order 10 approximation, and the one with short dashes is the order 20 approximation.Figure 2
Gram-Charlier approximations to the densityg0.3(x) (solid line) when q=1.4. The curve with long dashes is the order 10 approximation, and the one with short dashes is the order 20 approximation.Figure 3
Feasible region for the skewness and excess kurtosis of theGC(a,b;0,0,c3,c4) distributions and maximum likelihood estimate.


### 2.4. Exponential Change of Measure

The one-dimensional Cameron-Martin formula may be stated as follows: ifX~N(μ,σ2) then for q∈R and f≥0, (62)EeqXfX=eμq+1/2σ2q2EfX+σ2q.The same property may be expressed in terms of a change of measure. If X~PN(μ,σ2) and a change of measure is defined by(63)P′=eqXEPeqX·P,then X~P′N(μ+σ2q,σ2). The next result is an extension of the Cameron-Martin formula for the normal distribution to the Gram-Charlier distributions. There are different though equivalent formulas in Schlögl [10] (for the special case c1=c2=0).Theorem 4.
Suppose thatX~PGC(a,b;c1,…,cN) and that P′ is defined by (63) for q∈R. Then X~P′GC(a+b2q,b;c1′,…,cN′), where c1′,…,cN′ are found from (64)∑k=0Nbkck′sk=∑k=0Nbkckq+sk∑j=0Nbjcjqj,s∈R,or, more precisely, (65)ck′=1∑j=0Nbjcjqj∑l=kNlkbl-kclql-k.Proof.
It is sufficient to calculate the mgf ofX under P′: (66)EP′esX=1EPeqXEPeq+sX=⋯=ea+b2qs+1/2b2s2∑k=0Nbksk1∑j=0Nbjcjqj∑l=kNlkbl-kclql-k.Example 5.
IfX~PGC(a,b;0,0,c3,c4) and q=1 the change of measure (63) leads to (67)c1′=3b2c3+4b3c41+b3c3+b4c4,c2′=3bc3+6b2c41+b3c3+b4c4,c3′=c3+4bc41+b3c3+b4c4,c4′=c41+b3c3+b4c4.Hence, assuming c1=c2=0 does not imply that c1′ and c2′ are also zero. In other words, the family GC(a,b;0,0,c3,c4) is not closed under a change of measure that occurs naturally in option pricing (see Section 5).


### 2.5. Simulating Gram-Charlier Distributions

Simulation is required for many kinds of options, and it turns out that the Gram-Charlier distributions are very easy to generate, as we now show; there is no need to invert their distribution functions. When estimating some quantitym=Eg(X1,…,Xs) by simulation, one generates n independent vectors (68)X1j,…,Xsj,j=1,…,n,with the same distribution. Suppose all the X’s are independent and have the same Gram-Charlier distribution with density (69)fXx=ϕa,b,xpx,where the polynomial p(x) is given in (5). Then (70)EgX1,…,Xs=∫Rsgx1,…,xs∏k=1sϕa,b,xkpxkdx1⋯dxs=EQgX1,…,Xs∏k=1spXk,where under the measure Q the X’s have a normal distribution N(a,b2). Hence, estimating Eg(X1,…,Xs) by simulation can be performed by generating ordinary normal random variables X1j,…,Xsj,j=1,…,n, and then using the estimator (71)m^n=1n∑j=1ngX1j,…,Xsj∏k=1spXkj.This is an application of the likelihood ratio method.


## 2.1. TheGC(a,b;0,0,c3,c4) Family

Here the exact region for the(c3,c4) that lead to a true probability distribution has been found. This goes back to Barton and Dennis [16], but a more detailed explanation is given in Jondeau and Rockinger [6]. The region is shown in Figure 3 (use the correspondence s¯=6c3, k¯=24c4 from Theorem 2(h)). An important fact about this region is that it is not rectangular; the possible excess kurtosis values depend on skewness, and conversely.


## 2.2.GC(a,b;c~) Distributions of Order 4 and Higher

TheGC(a,b;c1,c2,c3,c4) family has six parameters, rather than four, and thus has more degrees of freedom in fitting data; to the authors’ knowledge the general six-parameter Gram-Charlier distribution has been used in financial applications by León et al. [11] only (those authors use the subfamily consisting of polynomials p(·) that are squares of some second-degree polynomial). Schlögl [10] fits the six-and eight-parameter families GC(a,b;0,0,c3,c4,c5,c6) and GC(a,b;0,0,c3,c4,c5,c6,c7,c8) to data.The set of(c1,c2,c3,c4) that yield true probability distributions has not been identified, but it is possible to fit the six parameters a,b,c1,c2,c3,c4 by maximum likelihood.


## 2.3. Why We Do Not Assume Thatc1=c2=0

Almost all previous authors have assumed thatc1=c2=0, because they used normalized data:(46)Y=X-EXVarX(see part (j) of Theorem 2) as explained below. In this section we explain why it is important not to restrict Gram-Charlier distributions or series in that way. The first reason for not doing so is that enlarging the parameter space can only be a good thing. The second one is that in fitting those distributions to data there may be very real advantages in letting c1 and c2 be different from zero. The only downside of letting c1 and c2 take nonzero values, and it is of no real importance, is that c3 and c4 lose their simple relationship with skewness and excess kurtosis (see part (h) of Theorem 2). A third reason is that after an exponential change of measure a Gram-Charlier distribution will rarely have c1=c2=0; see Section 2.4.We now show, using an example that can be worked out explicitly, that it is not always best to use normalized data when fitting Gram-Charlier distributions, because choosing another affine transformation of the data may well yield a much better fit.The “standard” Gram-Charlier expansion for a functionf(·) is(47)ϕx∑k=0∞ckHek,where (48)ck=1k!∫-∞∞Hekxfxdx=1k!EHekX.A classical result about Hermite series, proved by Cramér [17], is that sufficient conditions for the Gram-Charlier expansion (47) to converge to (49)12fx-+fx+for all x∈R are that (i) f(·) has finite variation in every bounded interval, and (ii) satisfies (50)∫-∞∞ex2/4fxdx<∞.If X has density f(·), then the last condition is EeX2/4<∞. This condition cannot be improved upon, in the sense that there are cases where the Gram-Charlier series defined above diverges, although EeρX2<∞ for all ρ<1/4 (it will be shown below that one such case is the normal distribution with mean 0 and variance 2).Let us first calculate(51)hq,k=EHekX,X~N0,q2.Using the generating function (16) we find (52)hq,k=0ifkis  oddk!2k/2k/2!q2-1k/2ifkis  even.The Gram-Charlier expansion corresponding to the N(0,q2) density is thus(53)ϕx∑k=0∞hq,kk!Hekx=ϕx∑n=0∞12nn!q2-1nHe2nx.It is possible to determine whether this converges or not when x=0; the Gram-Charlier series for f(0) is (54)12π∑n=0∞-1n2n!22nn!2q2-1n.From Stirling’s formula, (55)2n!22nn!2~1πnas n tends to infinity, which implies that there are three possibilities: expansion (54) (i) converges absolutely if 0<q2<2, (ii) converges simply if q2=2, and diverges if q2>2. (The preceding calculations are from Cramér [17].)Let us now consider a random variableXp, 0≤p≤1, with density (56)gpx=pϕx+1-p1pϕxp.In words, Xp has a N(0,1) distribution with probability p and a N(0,p2) distribution with probability 1-p. Define the normalized variable (57)Yp=XpVarXp=Xpσp,σp2=VarXp=p+p2-p3.The coefficients cp,k of the Gram-Charlier expansion for Yp are (58)cp,k=1k!EHekYp=1k!ph1/σp,k+1-php/σp,k.We have shown that (53) diverges if q2>2. Here we see that cp,k involves hq,k with q=1/σp and q=p/σp; in the second case there is no problem, as (p/σp)2<2 for all p∈(0,1]. In the first case, (59)1σp2>2⟺p+p2-p3<12.This condition is satisfied for p smaller than p∗≈0.403, and so the Gram-Charlier expansion for Yp diverges for all p<p∗. This is a sad state of affair: a series designed to work for distributions that are “close to the normal” that fails for combinations of two normal densities!There is an easy solution: use a different scaling forXp. Rather than multiplying Xp by 1/σp, use a factor r such that the series converges. In this example we know that the series converges if and only if r<2. Say we choose p=.3; this implies (60)1VarXp=10.3+0.32-0.33=2.75482>2.Figure 1 shows how the true density for X0.3 is approximated by order 10 and 20 Gram-Charlier truncated series if r=1/σ0.3≐1.66 is chosen; the graph is not surprising, the infinite Gram-Charlier expansion diverges. Figure 2 shows the same except that r is set to 1.4. The latter is smaller than 2, so the infinite expansion converges. In this example the density to be approximated is symmetrical about 0, so c1 is always 0, while (61)q=1σ0.3⟹c2=0;q=1.4⟹c2=-0.14426.Figure 1
Gram-Charlier approximations to the densityg0.3(x) (solid line) when q=1/σ0.3≐1.66, resulting in c1=c2=0. The curve with long dashes is the order 10 approximation, and the one with short dashes is the order 20 approximation.Figure 2
Gram-Charlier approximations to the densityg0.3(x) (solid line) when q=1.4. The curve with long dashes is the order 10 approximation, and the one with short dashes is the order 20 approximation.Figure 3
Feasible region for the skewness and excess kurtosis of theGC(a,b;0,0,c3,c4) distributions and maximum likelihood estimate.


## 2.4. Exponential Change of Measure

The one-dimensional Cameron-Martin formula may be stated as follows: ifX~N(μ,σ2) then for q∈R and f≥0, (62)EeqXfX=eμq+1/2σ2q2EfX+σ2q.The same property may be expressed in terms of a change of measure. If X~PN(μ,σ2) and a change of measure is defined by(63)P′=eqXEPeqX·P,then X~P′N(μ+σ2q,σ2). The next result is an extension of the Cameron-Martin formula for the normal distribution to the Gram-Charlier distributions. There are different though equivalent formulas in Schlögl [10] (for the special case c1=c2=0).Theorem 4.
Suppose thatX~PGC(a,b;c1,…,cN) and that P′ is defined by (63) for q∈R. Then X~P′GC(a+b2q,b;c1′,…,cN′), where c1′,…,cN′ are found from (64)∑k=0Nbkck′sk=∑k=0Nbkckq+sk∑j=0Nbjcjqj,s∈R,or, more precisely, (65)ck′=1∑j=0Nbjcjqj∑l=kNlkbl-kclql-k.Proof.
It is sufficient to calculate the mgf ofX under P′: (66)EP′esX=1EPeqXEPeq+sX=⋯=ea+b2qs+1/2b2s2∑k=0Nbksk1∑j=0Nbjcjqj∑l=kNlkbl-kclql-k.Example 5.
IfX~PGC(a,b;0,0,c3,c4) and q=1 the change of measure (63) leads to (67)c1′=3b2c3+4b3c41+b3c3+b4c4,c2′=3bc3+6b2c41+b3c3+b4c4,c3′=c3+4bc41+b3c3+b4c4,c4′=c41+b3c3+b4c4.Hence, assuming c1=c2=0 does not imply that c1′ and c2′ are also zero. In other words, the family GC(a,b;0,0,c3,c4) is not closed under a change of measure that occurs naturally in option pricing (see Section 5).


## 2.5. Simulating Gram-Charlier Distributions

Simulation is required for many kinds of options, and it turns out that the Gram-Charlier distributions are very easy to generate, as we now show; there is no need to invert their distribution functions. When estimating some quantitym=Eg(X1,…,Xs) by simulation, one generates n independent vectors (68)X1j,…,Xsj,j=1,…,n,with the same distribution. Suppose all the X’s are independent and have the same Gram-Charlier distribution with density (69)fXx=ϕa,b,xpx,where the polynomial p(x) is given in (5). Then (70)EgX1,…,Xs=∫Rsgx1,…,xs∏k=1sϕa,b,xkpxkdx1⋯dxs=EQgX1,…,Xs∏k=1spXk,where under the measure Q the X’s have a normal distribution N(a,b2). Hence, estimating Eg(X1,…,Xs) by simulation can be performed by generating ordinary normal random variables X1j,…,Xsj,j=1,…,n, and then using the estimator (71)m^n=1n∑j=1ngX1j,…,Xsj∏k=1spXkj.This is an application of the likelihood ratio method.


## 3. Convolution of Gram-Charlier Distributions; Gram-Charlier Processes

The simplest way to find the distribution of the sum of two independent Gram-Charlier distributed variables is to multiply their moment generating functions. SupposeXj~GC(aj,bj;c1(j),…,cN(j)), j=1,2 are independent. Then (72)EezX+Y=ea1+a2t+b12+b22t2/2∑k=0N1ck1tk∑k=0N2ck2tk.Expanding the product, this says that(73)X+Y~GCa1+a2,b12+b22;c1,…,cN1+N2,ck=∑j=0kcj1ck-j2.It is then possible to know the explicit distribution of (74)Zn=X1+⋯+Xn,where the X’s are independent and have a Gram-Charlier distribution, constituting a discrete-time Gram-Charlier process. If the X’s have the same GC(a,b;c1,…,cN) distribution then {Zn,n≥0} is a random walk. The derivation of the distribution of Zn can be done recursively, using (73), or it can be done by finding the Taylor expansion of (75)∑k=0Ncktkn=∑k=0nNckntk=1+c1nt+c2n+12c12nn-1t2+⋯+cNntnN,and thus (76)Zn~GCan,bn;c1n,c2n+12c12nn-1,…,cNn.These computations are simple using symbolic mathematics software. An example is given at the end of Section 6.1.The above raises the question of whether there is a continuous-time process that has Gram-Charlier distributed increments. There is such a process with normal increments (Brownian motion), and it is moreover a Lévy process.Theorem 6.
The only Lévy process with Gram-Charlier distributed increments is Brownian motion.Proof.
It is sufficient to show that, besides the normal distribution, any Gram-Charlier distribution cannot be infinitely divisible. IfX has a Gram-Charlier distribution then its moment generating function is of the form (77)Mz=EezX=eaz+b2z2/2πz,where π(·) is a polynomial. Suppose Xn,1,…,Xn,n are independent, have the same distribution, and add up to X (in law). Fubini’s theorem implies that EetXn,1 is finite for all real t, and thus (78)Mnz=EezXn,1is an analytic function of z in the whole complex plane. We then have the identity (79)Mnzn=Mz,z∈C,for every n=1,2,3,…. This means that the function π(z)1/n, which is well defined for all z that is not a zero of π, has an analytic continuation over the whole of C for each n≥1, which is impossible, for instance, take n larger than the degree of π.More precisely, this says that if we exclude Brownian motion, no increment of any Lévy process can have a Gram-Charlier distribution. Any Gram-Charlier process with independent increments must be discrete-time.


## 4. The Log Gram-Charlier Distribution

The distribution of the exponential of a Gram-Charlier distributed variable will naturally be called “log Gram-Charlier”, as we do for the lognormal: ifY~GC(a,b;c~), then L=exp⁡Y~LogGC(a,b;c~). The density of L is (80)fLz=1zfYlog⁡z,z>0.This distribution has all moments finite, and they are given by Theorem 2(b). The log Gram-Charlier distribution shares one property with the lognormal, it is “moment indeterminate.”Theorem 7.
The log Gram-Charlier distribution is not determined by its moments. More precisely, there is a noncountable number of other distributions that have the same moments as any particular log Gram-Charlier distribution.Proof.
There is a well-known way to construct a family of distributions that have the same moments as the lognormal (Feller [18], p. 227); the trick works for arbitrary parameters but, for simplicity, let L~LogN(0,1); that is, (81)fLz=1z2πe-1/2log2z,z>0.All the functions (82)fαz=fLz1+αsin⁡2πlog⁡z,-1≤α≤1,are nonnegative, integrate to 1, and have the same moments as L (as a result of the symmetry of the standard normal distribution, after an obvious change of variable). Now, suppose L~LogGC0,1;c~ (once again the same arguments work for other values of a and b), and consider (83)gαz=fLz+αϕzsin⁡2πlog⁡z.By the same change of variables used for the lognormal one immediately finds that gα integrates to one and has the same moments as L. The only difference with the case of the lognormal is that ga is not necessarily nonnegative for -1≤α≤1. Two cases may arise. The first one is that the polynomial (84)py=∑k=0NckHekyhas no real zero. In that case its infimum is strictly greater than zero, and one can find a nonempty interval I=(-ϵ,ϵ) such that gα is nonnegative for all α∈I. In the second case p(·) has one or more zeros and the previous argument breaks down but can be modified to yield the same conclusion, if the function sin⁡2πy is replaced with (85)sy=sin⁡2πy1y∉J,where J is a collection of intervals that include the zeros of p(·), so defined that s(·) is not identically zero and satisfies (86)sy+2kπ=sy,s-y=-sy. (These two conditions are sufficient for (87)∫-∞∞ekyϕysydy=0to hold for k=1,2,….) The details are omitted. Another way to prove that the log Gram-Charlier distribution is moment indeterminate when p(y) has no zero (and thus the density fL(z) does not take the value zero for any z>0) is to use a Krein condition (Stoyanov [19], p. 941), which says that a continuous distribution on R+ with positive density f(·) is not determined by its moments if(88)-∫0∞log⁡fz21+z2dz<∞.


## 5. Option Pricing Formulas

The formulas below hold for any vector(c1,…,cN) and thus extend those that have been derived by previous authors for the case where the log return has a GC(a,b;0,0,c3,c4) distribution or a “squared” Gram-Charlier distribution (León et al. [11]). Schlögl [10] derives a formula equivalent to (90) but for densities expressed as an infinite Gram-Charlier series with c1=c2=0.As previous authors have done, we consider a market with a risky securityS and a risk-free security with annual return r and suppose that the log return for period [0,T] (denoted as XT) has a Gram-Charlier distribution under the risk-neutral (or “pricing”) measure, which we denote as Q. The physical measure (usually denoted as P) is not specified (nothing says that the log return has or does not have a Gram-Charlier distribution under the physical measure).The market model may have one or more periods, but since we consider the pricing of ordinary European puts and calls only the distribution of the log return for the whole period[0,T] is required. For other types of options, in particular the applications presented in the next section, it may be necessary to use the one-period returns separately, as is done in Theorem 11.The risky security has priceS0 at time 0, and ST=S0eXT.Theorem 8.
Suppose that a risky security pays dividends at a constant rateδ over [0,T] and that the risk-free rate of interest is r. Suppose also that under the risk-neutral measure Q the log return of the risky security over [0,T] is XT~QGC(a,b;c1,…,cN), which satisfies the martingale condition(89)ea+b2/2∑k=0Nbkck=er-δT.Then the time-0 price of a European call option with maturity T and strike price K is(90)C0=S0e-δTΦd1-Ke-rTΦd2+Ke-rTϕd2∑k=1Nck∗Hek-1-d1-ckHek-1-d2,where d2, d1, and {ck∗} are given by (97), (101), and (105). The price of the corresponding European put is (91)P0=Ke-rTΦ-d2-S0e-δTΦ-d1-Ke-rTϕd2∑k=1Nck∗Hek-1-d1-ckHek-1-d2.If N=4, then the above simplify to (92)C0=S0e-δTΦd1-Ke-rTΦd2+bKe-rTϕd2c2+b-d2c3+b2-bd2+d22-1c4,P0=Ke-rTΦ-d2-S0e-δTΦ-d1-bKe-rTϕd2c2+b-d2c3+b2-bd2+d22-1c4.Proof.
Absence of arbitrage implies thatEQST=e(r-δ)TS0, or (93)S0ea+b2/2∑k=0Nbkck=er-δTS0.This justifies (89). The price of the call is (94)C0=e-rTEQST-K+=C0+-C0-C0+=e-rTEQST1ST>K,C0-=Ke-rTEQ1ST>K.The second part is easier to deal with (95)C0-=Ke-rTQXT>log⁡KS0.The probability of the event XT>log⁡(K/S0) can be calculated explicitly by recalling that (XT-a)/b~QGC(0,1;c1,…,cN) and using Theorem 3: (96)QXT>logKS0=QXT-ab>1blogKS0-a=Φd2+ϕd2∑k=1NckHek-1-d2,where(97)d2=1blog⁡S0K+a.To calculate C0+ use the exponential change of measure formula, defining (98)Q′=STEQST·Q=eXTEQeXT·Q.Then (99)EQST1ST>K=EQSTEQ′1ST>K.Since XT~Q′GC(a+b2,b;c1′,…,cN′),(100)EQSTEQ′1ST>K=S0er-δTQ′XT-a-b2b>1blog⁡KS0-a-b2=S0er-δTΦd1+S0er-δTϕd1∑k=1Nck′Hek-1-d1,where(101)ck′=1∑j=0Nbjcj∑l=kNlkbl-kcl,k=1,…,N,d1=1blogS0K+a+b2=d2+b.Hence, (102)C0=S0e-δTΦd1-Ke-rTΦd2+S0e-δTϕd1∑k=1Nck′Hek-1-d1-Ke-rTϕd2∑k=1NckHek-1-d2.Using(103)S0e-δTϕd1∑k=0Nbkck=S0ea+b2/2-rTϕd1=Ke-rTϕd2(a consequence of (89) and (101)), it is possible to write (104)S0e-δTϕd1ck′=Ke-rTϕd2ck∗,where(105)ck∗=ck′∑k=0Nbkck=∑l=kNlkbl-kcl,k=1,…,N.This proves (90).
The price of a European put can be found from the put-call parity identity(106)P0=C0-S0+Ke-rT.Finally, if N=4 then the summation in (90) becomes(107)∑k=14ck∗Hek-1-d1-∑k=14ckHek-1-d2=bc2+b2-bd2c3+b3-b2d2+bd22-bc4.This ends the proof.The option price formulas are of the form “Black-Scholes plus correction term.” Observe, however, that the values ofd1 and d2 are different from what they would be in the Black-Sholes formula. More precisely, in the Black-Scholes model (where ck=0 for all k≥1) we have a=EQXT and b2=VarQXT, but this does not happen with Gram-Charlier distributed log returns, first because of the martingale condition (89) and second because of the result in part (j) of Theorem 2.


### 5.1. Sensitivities (“Greeks”)

The previous literature includes formulas for sensitivities of Gram-Charlier option prices, but only in the case of the four-parameterGC(a,b;0,0,c3,c4) Gram-Charlier distributions; see Jurczenko et al. [7], Rouah and Vainberg [20], and Chateau [3]. Below we give sensitivities of the option prices calculated above with respect to all the parameters for a general Gram-Charlier distribution, taking the martingale condition (89) into account. This means that a is a function of δ,b,{ck},r and T (we might have written a=a(δ,b,{ck},r,T)). However, a is not a function of S0.Theorem 9.
LetC0 be the price of the European call option described in Theorem 8. Then (108)Δ=∂C0∂S0=e-δTΦd1+ϕd1∑k=1Nck′Hek-1-d1,γ=∂2C0∂S02=ea+b2/2-rTbS0ϕd1∑k=0NckHek-d2,ρ=∂C0∂r=KTe-rTΦd2+ϕd2∑k=1NckHek-1-d2,κ=∂C0∂b=S0e-δTϕd1∑k=0Nck′-c1′ck+1′+k+2ck+2′Hek-d1,∂C0∂cj=Ke-rTϕd2∑k=1j-1bkHej-1-k-d2-bj∑k=1Nck′Hek-1-d1,j≥1.In the formula for κ the symbols cN+1′ and cN+2′ are equal to zero.Proof.
The following lemma is obtained by elementary calculations.
Lemma  10.For any integrable random variableUand any constantK, (109)∂∂sEsU-K+=EU1sU>K.IfUhas a continuous densityfU, then(110)∂2∂s2EsU-K+=K2s3fUKs.
FromST=S0eXT and (100), (111)Δ=∂∂S0e-rTEQS0eXT-K+=1S0C0+=e-δTΦd1+ϕd1∑k=1Nck′Hek-1-d1.From the lemma, (112)γ=∂2∂2S0C0=e-rT∂2∂2S0EQS0U-K+=K2e-rTS03fUKS0,U=eX,where XT~QGC(a,b;c1,…,cN). Now the density of U may be expressed in terms of the density of G=(XT-a)/b~QGC(0,1;c1,…,cN): (113)fUu=1bufG1blog⁡u-a=1buϕy∑k=0NckHekyy=1/blog⁡u-a.Setting u=K/S0 and then using (103), we obtain (114)γ=Ke-rTbS02ϕd2∑k=0NckHek-d2=ea+b2/2-rTbS0ϕd1∑k=0NckHek-d2.Next, (115)ρ=∂∂rEQS0ea+bG-rT-Ke-rT+=∂∂rEQhG,r,where G~QGC(0,1,c~). For fixed G, the function g(r)=h(G,r) is absolutely continuous and thus (116)ρ=EQ∂∂rS0ea+bG-rT-Ke-rT1S0ea+bG-rT>Ke-rT.Now, from (89),(117)ea+bG-rT=ebG-b2/2-δT∑k=0Nbkckdoes not depend on r, so (118)ρ=EQKTe-rT1S0ea+bG-rT>Ke-rT=KTe-rTΦd2+ϕd2∑k=1NckHek-1-d2.The same line of reasoning shows that the sensitivity to b is (119)κ=S0e-rTEQ∂∂bea+bG1S0ea+bG>K.Define (120)mb=EebG-b2/2=∑k=0bkck,m′b=∑k=1kbk-1ck=c1∗.In order to shorten the formulas we will write (from (89)) (121)ea=er-δT-b2/2mband thus find an expression for (122)S0e-δTEQ∂∂bebG-b2/2mb1G>-d2=S0e-δT-b2/2mb2EQGebGmb-ebGm′b+bmb1G>-d2.For arbitrary y∈R, (123)EQebG1G>-y=eb2/2mbΦy+b+ϕy+b∑k=1Nck∗Hek-1-y-b.To calculate (124)EQGebG1G>-y=∂∂bEQebG1G>-ywe need (see (105)) (125)∂ck∗∂b=∂∂b∑l=kNlkbl-kcl=k+1ck+1∗,k=1,…,N-1,∂cN∗∂b=0.Then, noting that m′(b)=c1∗=c1′m(b) and ϕHek-1′=-ϕHek, (126)EQGebG1G>-y=bEQebG1G>-y+eb2/2mbc1′Φy+b+ϕy+b∑k=0N-2k+2ck+2′Hek-y-b+ϕy+b∑k=0Nck′Hek-y-b.Setting y=d2 and putting the pieces together yield the result. When c~=0 this reduces to S0e-δTϕ(d1), which is the sensitivity of the Black-Scholes call price to b=σT.
Finally, turn to the sensitivities with respect tocj, j≥1. First, let us derive ∂a/∂cj; from (89), (127)∂a∂cjea+b2/2∑k=0Nbkck+bjea+b2/2=0or∂a∂cj=-bj∑k=0Nbkck=-bjea+b2/2+δ-rT.Next,(128)e-rT∂∂cjEQS0ea+bG-K+=e-rT∂∂cj∫-d2∞S0ea+by-Kϕy∑k=0NckHekydy=e-rT∫-d2∞S0∂a∂cjea+byϕy∑k=0NckHekydy+e-rT∫-d2∞S0ea+by-KϕyHejydy.The first integral reduces to S0Δ∂a/∂cj (see the derivation of Δ above). The second integral is evaluated by repeated partial integration: if g(x) is continuously differentiable k times and does not grow too fast, (129)∫y∞gxϕxHejxdx=ϕy∑k=0j-1gkyHej-1-ky+∫y∞gjxϕxdx.Hence the second integral in (128) becomes (130)S0ea-bd2-rTϕd2∑k=1j-1bkHej-1-k-d2+S0ea-rT∫-d2∞bjebyϕydy=S0ea-bd2-rTϕd2∑k=1j-1bkHej-1-k-d2+bjS0ea-rT+b2/2Φd1.Putting all this together (using (103)) we get the formula for ∂C0/∂cj.


## 5.1. Sensitivities (“Greeks”)

The previous literature includes formulas for sensitivities of Gram-Charlier option prices, but only in the case of the four-parameterGC(a,b;0,0,c3,c4) Gram-Charlier distributions; see Jurczenko et al. [7], Rouah and Vainberg [20], and Chateau [3]. Below we give sensitivities of the option prices calculated above with respect to all the parameters for a general Gram-Charlier distribution, taking the martingale condition (89) into account. This means that a is a function of δ,b,{ck},r and T (we might have written a=a(δ,b,{ck},r,T)). However, a is not a function of S0.Theorem 9.
LetC0 be the price of the European call option described in Theorem 8. Then (108)Δ=∂C0∂S0=e-δTΦd1+ϕd1∑k=1Nck′Hek-1-d1,γ=∂2C0∂S02=ea+b2/2-rTbS0ϕd1∑k=0NckHek-d2,ρ=∂C0∂r=KTe-rTΦd2+ϕd2∑k=1NckHek-1-d2,κ=∂C0∂b=S0e-δTϕd1∑k=0Nck′-c1′ck+1′+k+2ck+2′Hek-d1,∂C0∂cj=Ke-rTϕd2∑k=1j-1bkHej-1-k-d2-bj∑k=1Nck′Hek-1-d1,j≥1.In the formula for κ the symbols cN+1′ and cN+2′ are equal to zero.Proof.
The following lemma is obtained by elementary calculations.
Lemma  10.For any integrable random variableUand any constantK, (109)∂∂sEsU-K+=EU1sU>K.IfUhas a continuous densityfU, then(110)∂2∂s2EsU-K+=K2s3fUKs.
FromST=S0eXT and (100), (111)Δ=∂∂S0e-rTEQS0eXT-K+=1S0C0+=e-δTΦd1+ϕd1∑k=1Nck′Hek-1-d1.From the lemma, (112)γ=∂2∂2S0C0=e-rT∂2∂2S0EQS0U-K+=K2e-rTS03fUKS0,U=eX,where XT~QGC(a,b;c1,…,cN). Now the density of U may be expressed in terms of the density of G=(XT-a)/b~QGC(0,1;c1,…,cN): (113)fUu=1bufG1blog⁡u-a=1buϕy∑k=0NckHekyy=1/blog⁡u-a.Setting u=K/S0 and then using (103), we obtain (114)γ=Ke-rTbS02ϕd2∑k=0NckHek-d2=ea+b2/2-rTbS0ϕd1∑k=0NckHek-d2.Next, (115)ρ=∂∂rEQS0ea+bG-rT-Ke-rT+=∂∂rEQhG,r,where G~QGC(0,1,c~). For fixed G, the function g(r)=h(G,r) is absolutely continuous and thus (116)ρ=EQ∂∂rS0ea+bG-rT-Ke-rT1S0ea+bG-rT>Ke-rT.Now, from (89),(117)ea+bG-rT=ebG-b2/2-δT∑k=0Nbkckdoes not depend on r, so (118)ρ=EQKTe-rT1S0ea+bG-rT>Ke-rT=KTe-rTΦd2+ϕd2∑k=1NckHek-1-d2.The same line of reasoning shows that the sensitivity to b is (119)κ=S0e-rTEQ∂∂bea+bG1S0ea+bG>K.Define (120)mb=EebG-b2/2=∑k=0bkck,m′b=∑k=1kbk-1ck=c1∗.In order to shorten the formulas we will write (from (89)) (121)ea=er-δT-b2/2mband thus find an expression for (122)S0e-δTEQ∂∂bebG-b2/2mb1G>-d2=S0e-δT-b2/2mb2EQGebGmb-ebGm′b+bmb1G>-d2.For arbitrary y∈R, (123)EQebG1G>-y=eb2/2mbΦy+b+ϕy+b∑k=1Nck∗Hek-1-y-b.To calculate (124)EQGebG1G>-y=∂∂bEQebG1G>-ywe need (see (105)) (125)∂ck∗∂b=∂∂b∑l=kNlkbl-kcl=k+1ck+1∗,k=1,…,N-1,∂cN∗∂b=0.Then, noting that m′(b)=c1∗=c1′m(b) and ϕHek-1′=-ϕHek, (126)EQGebG1G>-y=bEQebG1G>-y+eb2/2mbc1′Φy+b+ϕy+b∑k=0N-2k+2ck+2′Hek-y-b+ϕy+b∑k=0Nck′Hek-y-b.Setting y=d2 and putting the pieces together yield the result. When c~=0 this reduces to S0e-δTϕ(d1), which is the sensitivity of the Black-Scholes call price to b=σT.
Finally, turn to the sensitivities with respect tocj, j≥1. First, let us derive ∂a/∂cj; from (89), (127)∂a∂cjea+b2/2∑k=0Nbkck+bjea+b2/2=0or∂a∂cj=-bj∑k=0Nbkck=-bjea+b2/2+δ-rT.Next,(128)e-rT∂∂cjEQS0ea+bG-K+=e-rT∂∂cj∫-d2∞S0ea+by-Kϕy∑k=0NckHekydy=e-rT∫-d2∞S0∂a∂cjea+byϕy∑k=0NckHekydy+e-rT∫-d2∞S0ea+by-KϕyHejydy.The first integral reduces to S0Δ∂a/∂cj (see the derivation of Δ above). The second integral is evaluated by repeated partial integration: if g(x) is continuously differentiable k times and does not grow too fast, (129)∫y∞gxϕxHejxdx=ϕy∑k=0j-1gkyHej-1-ky+∫y∞gjxϕxdx.Hence the second integral in (128) becomes (130)S0ea-bd2-rTϕd2∑k=1j-1bkHej-1-k-d2+S0ea-rT∫-d2∞bjebyϕydy=S0ea-bd2-rTϕd2∑k=1j-1bkHej-1-k-d2+bjS0ea-rT+b2/2Φd1.Putting all this together (using (103)) we get the formula for ∂C0/∂cj.


## 6. Two Applications

### 6.1. Equity Indexed Annuities

Hardy [12] describes the main types of equity indexed annuities (EIAs), point-to-point, annual ratchet and high water mark, in which an embedded European call option has to be priced. We consider compound ratchet EIAs without life-of-contract guarantee to explore the dependence of ratchet premium options on skewness and kurtosis of returns.A single premiumP is paid by the policy-holder, and the benefit under the ratchet premium contract is ([12], p. 248) (131)Bn=P∏t=1n1+max⁡αStSt-1-1,0.Here, n is the term of the contract in years, α is the participation rate (a number between 0 and 1), and St is the value of the equity index, usually the S&P500 index.In [12], p. 249, the formula(132)Pe-r+αe-δΦd1-e-rΦd2nfor Bn is proved for the value of the premium option under a compound annual ratchet contract in the Black-Scholes model. The same reasoning will now be applied when index log returns have a Gram-Charlier distribution. To find a formula for the price of this option, we assume that under the risk-neutral measure the one-year log returns are independent Gram-Charlier distributed random variables (133)Rt~QGCa,b;c~,t=1,…,n.Independence of the variables across time implies that the value of the ratchet premium option is (if the annual risk-free rate of interest is r)(134)V0=EQe-rnBn=P∏t=1ne-rEQ1+max⁡αStSt-1-1,0.As Hardy notes ([12], p. 249), each factor in the product is the price of a one-year European call option on St, if initial index price and strike are both equal to one.Theorem 11.
The no-arbitrage price of the EIA ratchet premium option described above is(135)V0=Pe-r+αC0n,where C0 is the price of a one-year call (see (90)) with S0=1 and K=1. The first-order sensitivities are (136)∂V0∂ξ=αPe-r+αC0n-1∂C0∂ξ,with ξ replaced with one of S0, r, b, cj, and ∂C0/∂ξ given in Theorem 9. The second-order sensitivity with respect to S0 is (137)∂2C0∂S02=α2Pe-r+αC0n-2∂C0∂S02+αPe-r+αC0n-1∂2C0∂S02.We now show the effect of skewness and kurtosis on the ratchet premium option values; this is done using the four-parameter distributionGC(a,b;0,0,c3,c4), for which this is especially simple. We also compute prices using the six-parameter GC(a,b;c1,c2,c3,c4) distribution.Our application involves an annual ratchet, and so the optionC0 in Theorem 11 has a maturity of one year; we thus need the distribution of one-year returns. For pricing that is consistent with the (derivatives) market what is needed is the one-year distribution of returns under the risk-neutral measure, which should be obtained from observed derivative prices; see León et al. [11] or Rompolis and Tzavalis [21] for more details. For illustrative purposes, we estimated the four parameters of that distribution by maximum likelihood ([6, 22]) to obtain the parameters of the GC(a,b;0,0,c3,c4):(138)a=0.0687,b=0.1685,c3=-0.1150,c4=0.03598.This says that the skewness and excess kurtosis of the fitted distribution are (139)s¯=6c3=-0.6898,k¯=24c4=0.8634.This point is on the boundary of the feasible (s¯,k¯) region; see Figure 3.We use the parametersa and b while varying the values of c3,c4 to show how EIA premiums vary with skewness and excess kurtosis. Figure 4 shows the ratchet premium option values (Theorem 11) of an EIA with participation rate α=0.6 and single premium P=100, as a function of skewness s¯ and excess kurtosis k¯. The risk-free rate is 3%, the dividend rate is 2%, and the term is 7 years. The martingale condition (89) must hold and so a is replaced with a′ so that (89) is satisfied. (N.B. The Black-Scholes model corresponds to c3=c4=0, in which case a′=(r-δ)T-b2/2; this is -0.004193 with the values of r, δ, and b we are using and T=1.)Figure 4
Ratchet premium option prices as a function of skewness and excess kurtosis, withGC(a,b;0,0,c3,c4) distributions.The effect of varying(c3,c4) is quite significant. The highest and lowest ratchet premium options in the graph are 109.30 to 104.24, which is the range [95%,100%] as a proportion of the Black-Scholes premium option ($109.26), which has skewness and excess kurtosis equal to zero. (The maximum value of the premium option is reached at (c3,c4)=(0.0230,0.00332567) and the minimum at (-0.0836,0.163); these points correspond to (s¯,k¯) equal to (.138,.0798) and (0.023,0.003326), resp.) The dependence on (s¯,k¯) is nearly linear. This is easily explained. The price of a call is (140)C0=e-rTEQST-K=e-rT∫1/blog⁡K/S0-a′∞S0ea′+bx-Kϕx∑k=0NckHekxdx.The sensitivity of C0 to cj (j=3,4) is thus (141)e-rT∂a′∂cjEQST1ST>K+e-rT∫1/blog⁡K/S0-a′∞S0ea′+bx-KϕxHejxdx.The first term above is relatively small, and the second one is almost constant, since a′ does not change much with cj. This carries over to the sensitivity of the ratchet premium option V0 to changes in cj (see the first formula in Theorem 11).The parameters fitted by maximum likelihood for theGC(a,b;c1,c2,c3,c4) distribution are (142)a=0.1174,b=0.1595,c1=-0.3053675695201066,c2=0.09542079373489153,c3=-0.12383971126335243,c4=0.06120331530131559.The martingale condition implies a′=0.0451, and the EIA premium option has value 107.90. The skewness and excess kurtosis of the fitted distribution are s¯=-0.5437 and k¯=0.5092. (N.B. The feasible region for (s¯,k¯) in the GC(a,b;c1,c2,c3,c4) case is naturally larger than the one for the GC(a,b;0,0,c3,c4); in this case (s¯,k¯) are just outside the region in Figure 3.) It was checked numerically that the parameters c1 to c4 lead to a nonnegative function.Table1 shows ratchet premium option values based on Gram-Charlier with four parameters that are lower than the Black-Scholes based ones. Rows A and F are Black-Scholes cases, row B is our estimation of the four-parameter Gram-Charlier (as explained above), row C is the four-parameter Gram-Charlier with the largest negative skewness possible (see Figure 1), row D is similar but with the largest positive skewness, and row E has maximum kurtosis but no skewness. Rows A to E have the same parameter b, the one we estimated from S&P500 for a four-parameter Gram-Charlier, while row F has the same b as the six-parameter Gram-Charlier we estimated. Skewness and kurtosis are shown, followed by the break-even (or “fair”) α (“B-E α” in the table) and then the ratchet premium option values for two specific choices for α. Percentage differences with the values in row A (resp., F) are also displayed (“ΔA%” in the table). It is notable that in all cases the effect of moving away from Black-Scholes is more pronounced for the break-even participation rates than for the option values themselves. The scenario with the greatest impact (compared to Black-Scholes) is the one zero skewness and maximal excess kurtosis (row E). Insurers are unlikely to accept a cost of contract (V0) greater than the received initial premium ($100); they may thus wish to adjust the participation rate α down to the break-even rate (to wit, see Lin and Tan [23]). The Black-Scholes break-even participation rate in row A is 41.9%, but it increases from 44.3% up to 49.3% in rows B to E. These are significant increases, ranging from 5.5% to 17.5%. Rows F and G are based on two quite different risk-neutral probability distributions, but the ratchet premium options are surprisingly similar. The parameters we estimated may not be those that a practitioner would use, but what this example says is that skewness and kurtosis have a very real importance in the pricing of equity indexed annuities and that the fair participation rate is, at least here, more sensitive to nonzero skewness and excess kurtosis than are option values. (See also [24] regarding the separate effects of skewness and excess kurtosis on option values.)Table 1
Comparison of ratchet option values computed with the Black-Scholes formula and Gram-Charlier distributions. Rows A to F are based on the four-parameterGC(a,b;0,0,c3,c4), while row G is based on a six-parameter GC(a,b;c1,c2,c3,c4).
a ′ b c 1 c 2 c 3 c 4 s ¯ k ¯ B-Eα Δ ( A )% V 0   ( α = 0.6 ) Δ ( A )% V 0   ( α = 0.419 ) Δ ( A )% A −0.0042 0.1685 0 0 0 0 0 0 0.419 0 109.26 0 100.00 0 B −0.0035 0.1685 0 0 −0.1150 0.0360 −0.6898 0.8634 0.443 5.5 107.60 −1.5 98.92 −1.1 C −0.0034 0.1685 0 0 −0.1749 0.1021 −1.0493 2.4508 0.478 13.9 105.42 −3.5 97.50 −2.5 D −0.0051 0.1685 0 0 0.1749 0.1021 1.0493 2.4508 0.446 6.3 107.39 −1.7 98.78 −1.2 E −0.0043 0.1685 0 0 0 0.1667 0 4.0000 0.493 17.5 104.59 −4.3 96.96 −3.0  a ′ b c 1 c 2 c 3 c 4 s ¯ k ¯ B-Eα Δ ( F )% V 0   ( α = 0.6 ) Δ ( F )% V 0   ( α = 0.441 ) Δ ( F )%  F 0.0027 0.1595 0 0 0 .0 0 0 0.441 0 107.69 0 100.00 0 G 0.0451 0.1595 −.3054 .09542 −0.12384 0.06120 −0.5437 0.5091 0.438 −0.7 107.90 0.2 100.14 0.14Let us compare this with the approach in [11]. That paper is about a subclass of the Gram-Charlier distributions described in this paper; the authors call their distributions “seminonparametric” and have density (6), where p(·) is the square of the polynomial q(·) If q(x) is of order m, then it is possible to express p(x)=q(x)2 as(143)px=∑j=02mδjHejx.However, most nonnegative polynomials cannot be expressed as squares of another polynomial, so the order-m SNP family is a strict subset of the order-2m Gram-Charlier family. The order 2m SNP family has 2m+2 free parameters (the coefficients of an order-2m polynomial plus the location and dispersion parameters a,b, subject to density integrating to one, which in effect removes one parameter); by comparison, the GC(a,b;c1,…,c2m) distributions have 2m+2 parameters, but they are not “free” because the resulting density must be nonnegative. It is then nota priori clear which of SNP(2m) or GC(a,b;c1,…,c2m) would do best. We looked at this problem with the S&P500 annual prices, comparing SNP(2) and GC(a,b;0,0,c3,c4), which both have four parameters. The likelihood function is written the same way in both cases; there are just different restrictions on the parameters. The results for GC(a,b;0,0,c3,c4) are given in (138), and those for the SNP(2) are (144)a=-0.2445,b=0.1940,c1=1.6140,c2=1.1733,c3=0.4253,c4=0.07320.This translates into the following moments for the SNP(2) fitted distribution: (145)mean0.0687standard  deviation0.1671skewness-0.6289excess  kurtosis2.5502.The excess kurtosis is strikingly far from the data’s 0.8903. The likelihood function evaluated at that estimate was lower than the likelihood evaluated at the parameters of the GC(a,b;0,0,c3,c4) distribution. This is a puzzling result that does not mean that GC(a,b;0,0,c3,c4) will always do better than SNP(2); it does however show that restricting the search to squares of polynomials does have consequences. León et al. [11] derive formulas for European option prices, taking the martingale restriction into account, based on expression (143); in this paper we do the same calculation for an arbitrary Gram-Charlier distribution. It is easy to see that the SNP class is not closed under convolution; that is, the distribution of the sum of two independent SNP variables is in general not an SNP distribution; the general formulas we derive are essential in order to add independent Gram-Charlier distributed variables (including those in the SNP subset), and therefore the general formulas are needed to define Gram-Charlier processes. León et al. [11] are correct in pointing out that there is an advantage in using p(x)=q(x)2 as far as the nonnegativity constraint is concerned, but we believe that it is no more difficult to estimate the parameters of a general Gram-Charlier distribution, including the nonnegativity constraint in the maximization procedure (rather than limiting the parameter space). A full comparison of the numerical and statistical advantages/disadvantages of our approach versus the one in [11] is an interesting avenue for further research.


### 6.2. Lookback Options

Consider a European lookback call option with fixed strike, with payoff(S¯-K)+, where S¯ is the maximum of the stock price between time 0 and time T. We let the time points used to determine the maximum be equally spaced. There is no closed form formula for the distribution of the maximum price nor for the price of such a discretely monitored lookback option, so simulation is used. We generated normal random variables and used the technique described in Section 2.5. The assumptions are maturity T=1, monthly monitoring, annual volatility 0.1685 (where this number comes from is explained in the next subsection). The interest rate is set to 0, initial stock price is $100, and strike price is also $100.Figure5 shows the option prices for values of skewness and excess kurtosis that can be achieved with the GC(a,b;0,0,c3,c4) family for monthly rates of return. Roughly speaking, prices decrease when excess kurtosis increases, and also as kurtosis decreases, though the dependence on c3 and c4 becomes less linear as k¯ increases. When c3=c4=0 (Gaussian case) the option price is 11.277. The price is a little higher ($11.286) when s¯=0.75 and k¯=1.0. The lowest price ($10.163) occurs when s¯=0 and k¯=4. Skewness and kurtosis thus have a significant effect on lookback option prices; the effect will of course vary depending on the parameters. (Observe that all those prices are quite far from the Black-Scholes continuous-monitoring lookback price, $14.17.)Figure 5
Lookback option prices as a function of skewness and excess kurtosis.


## 6.1. Equity Indexed Annuities

Hardy [12] describes the main types of equity indexed annuities (EIAs), point-to-point, annual ratchet and high water mark, in which an embedded European call option has to be priced. We consider compound ratchet EIAs without life-of-contract guarantee to explore the dependence of ratchet premium options on skewness and kurtosis of returns.A single premiumP is paid by the policy-holder, and the benefit under the ratchet premium contract is ([12], p. 248) (131)Bn=P∏t=1n1+max⁡αStSt-1-1,0.Here, n is the term of the contract in years, α is the participation rate (a number between 0 and 1), and St is the value of the equity index, usually the S&P500 index.In [12], p. 249, the formula(132)Pe-r+αe-δΦd1-e-rΦd2nfor Bn is proved for the value of the premium option under a compound annual ratchet contract in the Black-Scholes model. The same reasoning will now be applied when index log returns have a Gram-Charlier distribution. To find a formula for the price of this option, we assume that under the risk-neutral measure the one-year log returns are independent Gram-Charlier distributed random variables (133)Rt~QGCa,b;c~,t=1,…,n.Independence of the variables across time implies that the value of the ratchet premium option is (if the annual risk-free rate of interest is r)(134)V0=EQe-rnBn=P∏t=1ne-rEQ1+max⁡αStSt-1-1,0.As Hardy notes ([12], p. 249), each factor in the product is the price of a one-year European call option on St, if initial index price and strike are both equal to one.Theorem 11.
The no-arbitrage price of the EIA ratchet premium option described above is(135)V0=Pe-r+αC0n,where C0 is the price of a one-year call (see (90)) with S0=1 and K=1. The first-order sensitivities are (136)∂V0∂ξ=αPe-r+αC0n-1∂C0∂ξ,with ξ replaced with one of S0, r, b, cj, and ∂C0/∂ξ given in Theorem 9. The second-order sensitivity with respect to S0 is (137)∂2C0∂S02=α2Pe-r+αC0n-2∂C0∂S02+αPe-r+αC0n-1∂2C0∂S02.We now show the effect of skewness and kurtosis on the ratchet premium option values; this is done using the four-parameter distributionGC(a,b;0,0,c3,c4), for which this is especially simple. We also compute prices using the six-parameter GC(a,b;c1,c2,c3,c4) distribution.Our application involves an annual ratchet, and so the optionC0 in Theorem 11 has a maturity of one year; we thus need the distribution of one-year returns. For pricing that is consistent with the (derivatives) market what is needed is the one-year distribution of returns under the risk-neutral measure, which should be obtained from observed derivative prices; see León et al. [11] or Rompolis and Tzavalis [21] for more details. For illustrative purposes, we estimated the four parameters of that distribution by maximum likelihood ([6, 22]) to obtain the parameters of the GC(a,b;0,0,c3,c4):(138)a=0.0687,b=0.1685,c3=-0.1150,c4=0.03598.This says that the skewness and excess kurtosis of the fitted distribution are (139)s¯=6c3=-0.6898,k¯=24c4=0.8634.This point is on the boundary of the feasible (s¯,k¯) region; see Figure 3.We use the parametersa and b while varying the values of c3,c4 to show how EIA premiums vary with skewness and excess kurtosis. Figure 4 shows the ratchet premium option values (Theorem 11) of an EIA with participation rate α=0.6 and single premium P=100, as a function of skewness s¯ and excess kurtosis k¯. The risk-free rate is 3%, the dividend rate is 2%, and the term is 7 years. The martingale condition (89) must hold and so a is replaced with a′ so that (89) is satisfied. (N.B. The Black-Scholes model corresponds to c3=c4=0, in which case a′=(r-δ)T-b2/2; this is -0.004193 with the values of r, δ, and b we are using and T=1.)Figure 4
Ratchet premium option prices as a function of skewness and excess kurtosis, withGC(a,b;0,0,c3,c4) distributions.The effect of varying(c3,c4) is quite significant. The highest and lowest ratchet premium options in the graph are 109.30 to 104.24, which is the range [95%,100%] as a proportion of the Black-Scholes premium option ($109.26), which has skewness and excess kurtosis equal to zero. (The maximum value of the premium option is reached at (c3,c4)=(0.0230,0.00332567) and the minimum at (-0.0836,0.163); these points correspond to (s¯,k¯) equal to (.138,.0798) and (0.023,0.003326), resp.) The dependence on (s¯,k¯) is nearly linear. This is easily explained. The price of a call is (140)C0=e-rTEQST-K=e-rT∫1/blog⁡K/S0-a′∞S0ea′+bx-Kϕx∑k=0NckHekxdx.The sensitivity of C0 to cj (j=3,4) is thus (141)e-rT∂a′∂cjEQST1ST>K+e-rT∫1/blog⁡K/S0-a′∞S0ea′+bx-KϕxHejxdx.The first term above is relatively small, and the second one is almost constant, since a′ does not change much with cj. This carries over to the sensitivity of the ratchet premium option V0 to changes in cj (see the first formula in Theorem 11).The parameters fitted by maximum likelihood for theGC(a,b;c1,c2,c3,c4) distribution are (142)a=0.1174,b=0.1595,c1=-0.3053675695201066,c2=0.09542079373489153,c3=-0.12383971126335243,c4=0.06120331530131559.The martingale condition implies a′=0.0451, and the EIA premium option has value 107.90. The skewness and excess kurtosis of the fitted distribution are s¯=-0.5437 and k¯=0.5092. (N.B. The feasible region for (s¯,k¯) in the GC(a,b;c1,c2,c3,c4) case is naturally larger than the one for the GC(a,b;0,0,c3,c4); in this case (s¯,k¯) are just outside the region in Figure 3.) It was checked numerically that the parameters c1 to c4 lead to a nonnegative function.Table1 shows ratchet premium option values based on Gram-Charlier with four parameters that are lower than the Black-Scholes based ones. Rows A and F are Black-Scholes cases, row B is our estimation of the four-parameter Gram-Charlier (as explained above), row C is the four-parameter Gram-Charlier with the largest negative skewness possible (see Figure 1), row D is similar but with the largest positive skewness, and row E has maximum kurtosis but no skewness. Rows A to E have the same parameter b, the one we estimated from S&P500 for a four-parameter Gram-Charlier, while row F has the same b as the six-parameter Gram-Charlier we estimated. Skewness and kurtosis are shown, followed by the break-even (or “fair”) α (“B-E α” in the table) and then the ratchet premium option values for two specific choices for α. Percentage differences with the values in row A (resp., F) are also displayed (“ΔA%” in the table). It is notable that in all cases the effect of moving away from Black-Scholes is more pronounced for the break-even participation rates than for the option values themselves. The scenario with the greatest impact (compared to Black-Scholes) is the one zero skewness and maximal excess kurtosis (row E). Insurers are unlikely to accept a cost of contract (V0) greater than the received initial premium ($100); they may thus wish to adjust the participation rate α down to the break-even rate (to wit, see Lin and Tan [23]). The Black-Scholes break-even participation rate in row A is 41.9%, but it increases from 44.3% up to 49.3% in rows B to E. These are significant increases, ranging from 5.5% to 17.5%. Rows F and G are based on two quite different risk-neutral probability distributions, but the ratchet premium options are surprisingly similar. The parameters we estimated may not be those that a practitioner would use, but what this example says is that skewness and kurtosis have a very real importance in the pricing of equity indexed annuities and that the fair participation rate is, at least here, more sensitive to nonzero skewness and excess kurtosis than are option values. (See also [24] regarding the separate effects of skewness and excess kurtosis on option values.)Table 1
Comparison of ratchet option values computed with the Black-Scholes formula and Gram-Charlier distributions. Rows A to F are based on the four-parameterGC(a,b;0,0,c3,c4), while row G is based on a six-parameter GC(a,b;c1,c2,c3,c4).
a ′ b c 1 c 2 c 3 c 4 s ¯ k ¯ B-Eα Δ ( A )% V 0   ( α = 0.6 ) Δ ( A )% V 0   ( α = 0.419 ) Δ ( A )% A −0.0042 0.1685 0 0 0 0 0 0 0.419 0 109.26 0 100.00 0 B −0.0035 0.1685 0 0 −0.1150 0.0360 −0.6898 0.8634 0.443 5.5 107.60 −1.5 98.92 −1.1 C −0.0034 0.1685 0 0 −0.1749 0.1021 −1.0493 2.4508 0.478 13.9 105.42 −3.5 97.50 −2.5 D −0.0051 0.1685 0 0 0.1749 0.1021 1.0493 2.4508 0.446 6.3 107.39 −1.7 98.78 −1.2 E −0.0043 0.1685 0 0 0 0.1667 0 4.0000 0.493 17.5 104.59 −4.3 96.96 −3.0  a ′ b c 1 c 2 c 3 c 4 s ¯ k ¯ B-Eα Δ ( F )% V 0   ( α = 0.6 ) Δ ( F )% V 0   ( α = 0.441 ) Δ ( F )%  F 0.0027 0.1595 0 0 0 .0 0 0 0.441 0 107.69 0 100.00 0 G 0.0451 0.1595 −.3054 .09542 −0.12384 0.06120 −0.5437 0.5091 0.438 −0.7 107.90 0.2 100.14 0.14Let us compare this with the approach in [11]. That paper is about a subclass of the Gram-Charlier distributions described in this paper; the authors call their distributions “seminonparametric” and have density (6), where p(·) is the square of the polynomial q(·) If q(x) is of order m, then it is possible to express p(x)=q(x)2 as(143)px=∑j=02mδjHejx.However, most nonnegative polynomials cannot be expressed as squares of another polynomial, so the order-m SNP family is a strict subset of the order-2m Gram-Charlier family. The order 2m SNP family has 2m+2 free parameters (the coefficients of an order-2m polynomial plus the location and dispersion parameters a,b, subject to density integrating to one, which in effect removes one parameter); by comparison, the GC(a,b;c1,…,c2m) distributions have 2m+2 parameters, but they are not “free” because the resulting density must be nonnegative. It is then nota priori clear which of SNP(2m) or GC(a,b;c1,…,c2m) would do best. We looked at this problem with the S&P500 annual prices, comparing SNP(2) and GC(a,b;0,0,c3,c4), which both have four parameters. The likelihood function is written the same way in both cases; there are just different restrictions on the parameters. The results for GC(a,b;0,0,c3,c4) are given in (138), and those for the SNP(2) are (144)a=-0.2445,b=0.1940,c1=1.6140,c2=1.1733,c3=0.4253,c4=0.07320.This translates into the following moments for the SNP(2) fitted distribution: (145)mean0.0687standard  deviation0.1671skewness-0.6289excess  kurtosis2.5502.The excess kurtosis is strikingly far from the data’s 0.8903. The likelihood function evaluated at that estimate was lower than the likelihood evaluated at the parameters of the GC(a,b;0,0,c3,c4) distribution. This is a puzzling result that does not mean that GC(a,b;0,0,c3,c4) will always do better than SNP(2); it does however show that restricting the search to squares of polynomials does have consequences. León et al. [11] derive formulas for European option prices, taking the martingale restriction into account, based on expression (143); in this paper we do the same calculation for an arbitrary Gram-Charlier distribution. It is easy to see that the SNP class is not closed under convolution; that is, the distribution of the sum of two independent SNP variables is in general not an SNP distribution; the general formulas we derive are essential in order to add independent Gram-Charlier distributed variables (including those in the SNP subset), and therefore the general formulas are needed to define Gram-Charlier processes. León et al. [11] are correct in pointing out that there is an advantage in using p(x)=q(x)2 as far as the nonnegativity constraint is concerned, but we believe that it is no more difficult to estimate the parameters of a general Gram-Charlier distribution, including the nonnegativity constraint in the maximization procedure (rather than limiting the parameter space). A full comparison of the numerical and statistical advantages/disadvantages of our approach versus the one in [11] is an interesting avenue for further research.


## 6.2. Lookback Options

Consider a European lookback call option with fixed strike, with payoff(S¯-K)+, where S¯ is the maximum of the stock price between time 0 and time T. We let the time points used to determine the maximum be equally spaced. There is no closed form formula for the distribution of the maximum price nor for the price of such a discretely monitored lookback option, so simulation is used. We generated normal random variables and used the technique described in Section 2.5. The assumptions are maturity T=1, monthly monitoring, annual volatility 0.1685 (where this number comes from is explained in the next subsection). The interest rate is set to 0, initial stock price is $100, and strike price is also $100.Figure5 shows the option prices for values of skewness and excess kurtosis that can be achieved with the GC(a,b;0,0,c3,c4) family for monthly rates of return. Roughly speaking, prices decrease when excess kurtosis increases, and also as kurtosis decreases, though the dependence on c3 and c4 becomes less linear as k¯ increases. When c3=c4=0 (Gaussian case) the option price is 11.277. The price is a little higher ($11.286) when s¯=0.75 and k¯=1.0. The lowest price ($10.163) occurs when s¯=0 and k¯=4. Skewness and kurtosis thus have a significant effect on lookback option prices; the effect will of course vary depending on the parameters. (Observe that all those prices are quite far from the Black-Scholes continuous-monitoring lookback price, $14.17.)Figure 5
Lookback option prices as a function of skewness and excess kurtosis.


## 7. Discussion and Conclusion

This paper tries to provide a framework for applying Gram-Charlier series to option pricing. We have shown that option prices may be significantly affected by varying skewness or kurtosis away from their Gaussian values. The practical consequence is that using the Black-Scholes formula (i.e., zero skewness and excess kurtosis) may distort option prices. Gram-Charlier distributions capture skewness and kurtosis and retain a lot of the tractability of the normal distribution. Estimating parameters remains a challenge for Gram-Charlier distributions involving more than four parameters.Convergence of Gram-Charlier and other expansions in option pricing and other applied problems is an area for further research. The heuristic derivation of the Gram-Charlier/Edgeworth series of a densityg(·) in terms of the normal density and its derivatives (see, e.g., [25], Section  3.2) may give the impression that the inverted series should converge or at least have some asymptotic property; this is not the case in general. Gram-Charlier (or Edgeworth) expansions do not always converge to the true probability distribution, as Section 2.3 indicates; see [26], Section  17, or [25] for more details. In the context of the Central Limit Theorem, the Edgeworth series for the density of the normalized sum or random variables have asymptotic properties, when the number of summands n tends to infinity. Details are given in [25]. Multivariate versions of the Gram-Charlier series exist, for instance, [10, 27, 28], apply such generalizations.The Gram-Charlier type B series are based on the Poisson distribution, rather than the normal. There are several other expansions that can be used. For instance, series of Laguerre polynomials may be used for densities on the positive half-line (a convergent Laguerre series has some numerical success in pricing Asian options; see [29]). The Jacobi polynomials [30] are orthogonal over the interval (-1,1), but any other finite interval may be chosen by making the obvious change of variable, and this leads to series expansions for densities with compact support, based on the beta distribution. An application of log-transformed shifted Jacobi polynomials to arbitrary approximation of probability distributions on the positive half-line may be found in [31].



---

*Source: 8690491-2017-02-08.xml*