MODEL-FREE IMPLIED VOLATILITY: FROM SURFACE TO INDEX
M. FUKASAWA, I. ISHIDA, N. MAGHREBI, K. OYA, M. UBUKATA, AND K. YAMAZAKI
Abstract. We propose a new method for approximating the expected quadratic variation of anasset based on its option prices. The quadratic variation of an asset price is often regarded asa measure of its volatility, and its expected value under pricing measure can be understood as
the market’s expectation of future volatility. We utilize the relation between the asset varianceand the Black-Scholes implied volatility surface, and discuss the merits of this new model-free
approach compared to the CBOE procedure underlying the VIX index. The interpolation scheme
for the volatility surface we introduce is designed to be consistent with arbitrage bounds. We
show numerically under the Heston stochastic volatility model that this approach significantly
reduces the approximation errors, and we further provide empirical evidence from the Nikkei 225
options that the new implied volatility index is more accurate in predicting future volatility.

Key words: model-free implied volatility index; volatility forecasting; volatility surface; variance swaps
JEL classification: C02

1. Introduction
Market volatility is one of the most important ingredients in financial decision-making. Yet it is nota quantity directly observable in the market. One measure of market volatility is the options-implied
expectation of the short-term quadratic variation of the market index. The VIX index disseminated by
the Chicago Board Options Exchange (CBOE) is one such measure calculated from the prices of the S&P
500 index options. This volatility index is model-free in the sense that it does not rely on a particular option
pricing model such as the Black-Scholes model. Indeed, it approximates the risk-neutral expectationsof the quadratic variations of the market index over a fixed expiration period based on their model-freelink with the prices of the corresponding index options. The theory behind this link assumes, among
other things, the existence of a continuum of options over an infinite range of strike prices. Thus,
approximation errors would necessarily ensue in any attempt to quantify market volatility from thelimited number of observable option prices.
In order to reduce approximation errors, this paper proposes a new method based on a formula
induced by the model-free link introduced in the context of pricing variance swaps. Given a continuous
Electronic version of an article published as
International Journal of Theoretical and Applied Finance Vol. 14, No. 4 (2011) 433-463
DOI: 10.1142/S0219024911006681
c World Scientific Publishing Company
http://www.worldscinet.com/ijtaf/ .
1

2

M. FUKASAWA, I. ISHIDA, N. MAGHREBI, K. OYA, M. UBUKATA, AND K. YAMAZAKI

semimartingale S standing for an asset price process, we have
∫
1
(1.1)
E[hlog(S)iT ] =σ(g(z))2 φ(z)dz.
T
Here, σ expresses the Black-Scholes implied volatility of the asset S with maturity T as a function oflog-moneyness k = log(K/F) with strike price K and forward price F, g is the inverse function of the
mapping
(1.2)
where
(1.3)

k 7→ d2 (k) := d2 (k, σ(k))√
k
σ T
d2 (k, σ) := − √ −,2σ Tand φ is the standard normal density. It should be noted that the use of the Black-Scholes implied
volatility function σ does not imply the dependency of the new method on this particular option theory.
The Black-Scholes implied volatility function is a nonlinear transformation of the option prices that reduces
to a constant only in a special case where the underlying asset follows a log-normal distribution. The
formula (1.1) states that an average of the Black-Scholes implied variance σ2 is equal to an average of thequadratic variations. This issue will be discussed in further detail in Section 2.
The form of the integral (1.1) written with respect to the standard normal density φ, together withpolynomial interpolation of the integrand, enables us to avoid numerical integration. Our approximation
method therefore provides better accuracy than those which rely on discretization to compute integralssuch as the CBOE method and the one proposed by Jiang and Tian (2007). We also propose a new
C1 -interpolation method by cubic polynomials in the Black-Scholes implied volatility scale to overcome
the errors typically incurred by the usual C2 -interpolation in terms of cubic splines.
We develop a new volatility index using our approximation method and evaluate its eﬃciency numerically based on the Nikkei 225 index options traded on the Osaka Securities Exchange (OSE). Before
using the actual options price data, we first consider artificial data generated under the Heston stochastic volatility model (Heston (1993)). The result shows that the new index values are significantly closer
to the true values than those produced by the CBOE procedure. We then use the actual options dataand evaluate empirically its ability to forecast market volatility, using as many as 15 realized varianceestimators of the Nikkei 225 stock average. Our new index is associated with higher R2 ’s than its CBOE
counterpart for all cases, implying that it forecasts market volatility more accurately. These numericalresults provide further evidence on the eﬀectiveness of our approximation method and the usefulnessof the new volatility index.
The remainder of the paper is organized as follows. We first review in Section 2 the theory underlying
the development of the model-free implied volatility, and describe in Section 3 the new algorithm forcomputing the volatility index. We then discuss the advantages of using our approximation method
over the CBOE procedure in Section 4, along with numerical results within the setting of the Heston
model. Section 5 provides some empirical evidence based on the market prices of the Nikkei 225
options, and shows that the new volatility index has the potential to predict future volatility moreaccurately than its comparative index based on the CBOE procedure. Section 6 concludes the paper.

MODEL-FREE IMPLIED VOLATILITY: FROM SURFACE TO INDEX

3

For readers’ convenience, Appendix A gives the proofs of the theorems in Section 2 that are extracted
from Fukasawa (2010). Appendix B provides a detailed description of the various estimators of realized
variance used in Section 5.
2. Model-free formulas and their approximations
This section reviews the model-free formulas for the quadratic variation of the underlying asset andthe approximation methods. Let S and S0 stand for risky and risk-free asset price processes, respectively.
If S is a continuous semimartingale and S0 is a deterministic process of locally bounded variation, thenwe have, by Itô’s formula,
 0
∫
∫ T
( ) ∫ TdS∗t ST  1
ST
dSt 1 T dhSit  − hlog(S)iT
log
=
−
=
+
log
 0 2
∗
S0
St
2 0
St
S2t
S0
0
0
where S∗t := St /S0t . This implies, under any risk-neutral expectation E, that
 0 
[
(
)]
[
( )]

 ST 
[
]
ST
ST




(2.1)
E hlog(S)iT = −2E log+ 2E log  0  = −2E log.
S0
E[ST ]S0
The expected quadratic variation is therefore determined by the distribution of ST /E[ST ]. It should beremarked here that this relation between the expected quadratic variation and the expected log price
was extended to the time-changed Lévy model by Carr, Lee and Wu (2009), where the multiplier in theabove equality changes from −2 in general.
Now, denote by F the T-expiry forward price of the asset S, and by C(K) and P(K), respectively, theundiscounted call and put option prices of the same asset with strike K and maturity T. As usual, we
suppose that there exists a probability measure E such thatF = E[ST ], C(K) = E[(ST − K)+ ], P(K) = E[(K − ST )+ ]for all K > 0. In addition, we suppose that ST admits a density under E and that there exists p > 0 suchthat−p

E[ST ] < ∞.
Then, the density apparently coincides with C00 (K) = P00 (K), and after integration-by-parts we have for
any K0 > 0 that
∫ ∞
∫ F
∫ K0
C(K)
P(K)
K−F
dK + 2dK + 2dK.
(2.2)−2E[log(ST /F)] = 2
2
2
2
K
K
K0
K0 K
0
Combining this and (2.1), we obtain a model-free relation between the expected quadratic variations andthe vanilla option prices, which enables us to calculate a model-free measure of the volatility implicit in
the market prices of the options contracts. In practice, however, the number of available option prices
is finite and we need to approximate the integrals on the right-hand side of (2.2). The CBOE procedure,
which underlies the computation of the current VIX index, refers to the approximation formula:(2.3)−2E[log(ST /F)] ≈ 2

K∑
0 −1
K=Kmin

K
)2
(
max
∑
P(K)
P(K0 ) + C(K0 )
C(K)
F
−1 .
∆K +
∆K0 + 2
∆K −
K0
K2
K2
K02
K=K +10

4

M. FUKASAWA, I. ISHIDA, N. MAGHREBI, K. OYA, M. UBUKATA, AND K. YAMAZAKI

Here, Kmin and Kmax represent, respectively, the lowest and highest strikes among the listed options.
The strike price K0 is taken to be closest to the forward price F among the set of strike prices smaller
than or equal to F. This approximation method will be revisited in Section 4, in comparison to ourown approximation method described in the next section. For a comprehensive account of the CBOE
procedure, reference can be made to, for example, Carr and Wu (2006).
The formula (1.1) is given in Gatheral (2006) as a reformulation of (2.2). See also Carr and Lee (2009).
The essentially same formula can be found in Morokoﬀ, Akesson and Zhou (1999) and Chriss andMorokoﬀ (1999). Here we state a result from Fukasawa (2010) where the formula is proved under lessrestrictive conditions. We give the proof in Appendix A.
Theorem 2.1. Define PBS : R × (0, ∞) → (0, ∞) by
PBS (k, σ) := Fek Φ(−d2 (k, σ)) − FΦ(−d1 (k, σ)),√where Φ is the standard normal distribution function, d2 (·, ·) is as defined in (1.3), and d1 (k, σ) := d2 (k, σ) + σ T.
Then the Black-Scholes implied volatility σ : R → [0, ∞) is well-defined by
σ(k) := PBS (k, ·)−1 (P(Fek )),(2.4)
or, equivalently,

PBS (k, σ(k)) = P(Fek ).
Moreover, the mapping d2 : k 7→ d2 (k, σ(k)) is a decreasing function and it holds that∫ ∞2
− E[log(ST /F)] =σ(g(z))2 φ(z)dz,
T−∞
where g is the inverse function of d2 .
As we shall shortly see in the next section, we successfully avoid numerical integrations by taking
advantage of the form of the integral with respect to φ and the polynomial interpolation of the integrand.
As can be easily seen in the analysis in Appendix A, the condition that ST admits a density also implies
that σ is continuously diﬀerentiable and its derivative is absolutely continuous. It is therefore natural
to take a C1 -interpolation/extrapolation by a piecewise polynomial. In addition, we have the following
result that should be taken into consideration. For its proof, see Appendix A.
Theorem 2.2. For z, z0 ∈ R, put
σ̂(z) := σ(g(z)),√α± (z; z0 ) := −z ± σ̂(z0 )2 + 2z0 σ̂(z0 ) + z2 .
Then we have the following.
(i) For every z ∈ R, we havedσ̂(z) > −1.
dz(ii) For every z, z0 ∈ R with z > z0 ≥ 0, we haveσ̂(z) > α+ (z, z0 ) > (σ̂(z0 ) + z0 − z)+ .

MODEL-FREE IMPLIED VOLATILITY: FROM SURFACE TO INDEX

5(iii) For every z, z0 ∈ R with z0 > z ≥ 0, we haveσ̂(z) < α+ (z, z0 ) < σ̂(z0 ) + z0 − z.
(iv) For every z, z0 ∈ R with z < z0 ≤ 0, we haveσ̂(z) > α− (z, z0 ).
(v) There exists z∗ < 0 such that σ̂(z∗ ) = −z∗ and it holds for every z, z0 ∈ R with z < z0 ≤ z∗ that−z > σ̂(z) > α− (z, z0 ) > 0.
The principal message from this theorem is the following:(i) The interpolation scheme should not produce excessive oscillations.
(ii) The extrapolation scheme should not induce a rapid decay of σ̂(z) as |z| → ∞. In fact Theorem 2.2
implies that 1/σ̂(z) = O(|z|) as |z| → ∞.
In the next section, we propose an algorithm that meets the above requirements. We further address its
advantages over the CBOE procedure in Section 4.
3. Algorithm
This section describes our algorithm for approximating accurately the annualized expected quadratic
variation
∫
]
2
1 [
(3.1)
E hlog(S)iT = − E[log(ST /F)] =σ(g(z))2 φ(z)dz
TTusing a set of option prices of the same maturity. It is applied to several options data sets with diﬀerent
maturities, and the results are then combined to attain the volatility index. In a nutshell, the algorithm
consists of three steps: (1) selecting a set of valid options to be used, (2) approximating the function
z 7→ σ(g(z))2 in terms of a set of cubic polynomials and (3) integrating the obtained function with respect
to the normal density φ. The rest of this section is devoted to describing each step in detail.
3.1. STEP 1 – selecting the options to be used in the index calculation. As a preliminary step, a subset
of options are selected for the index calculation. We first identify the at-the-money (ATM) strike price K0 ,and subsequently the out-of-the-money (OTM) calls and puts. The forward price is determined throughthe theoretical put-call parity relationship using ATM put and call options. Both transaction and bid/ask
data are used in this step.
The first task of this step is to determine K0 , which corresponds to the strike price at which thediﬀerence between put and call option prices is minimized. Here, we use only transactions data. Notice,
however, that transactions prices are used only for the purpose of determining K0 , thereby contributing
only indirectly to the index calculation. In case of equal put-call diﬀerences between two sets of options,the highest strike is selected. The forward price is then calculated using the formula
F = K0 + erT (call price − put price)where the call and put option prices are ATM transaction prices with strike K0 and we use the 3-month
Certificates of Deposit (CD) rates as proxy for the interest rate r.

6

M. FUKASAWA, I. ISHIDA, N. MAGHREBI, K. OYA, M. UBUKATA, AND K. YAMAZAKI

Using K0 as the ATM strike price, the sets of OTM puts and calls are then identified. OTM puts are
selected among all put options with strike prices less than or equal to K0 , whereas OTM calls are selectedamong all call options with strike prices greater than K0 . We consider only options for which bid/ask
prices are well-defined. This means that an option is bound to be eliminated unless both bid and ask
prices are available. We further sift the data by discarding any option with a ratio of ask to bid prices
greater than or equal to some constant c so as to focus on the set of reliable mid quotes that will be usedlater. Here we use c = 2. It should be mentioned that the options at K0 themselves may fail to satisfy
these conditions and may not be used for the index calculation.
Table 1 illustrates how the ATM strike price and OTM puts and calls are selected. The strike price
10, 000 (in JPY) is associated with the smallest diﬀerence 105 (using transaction data) and it hence becomes
the ATM strike K0 . The puts and calls are then selected accordingly using the bid-ask data. Here theput options with strike 6, 500 and 7, 500 and the call options with 12, 500 and 12, 750 are eliminated by
the ratio rule described above. Suppose r = 0.004825, which represents the interest rate over the time tomaturity T = 0.11984398782344, expressed in years, then the forward price F amounts to10000 + exp(0.004825 × 0.11984398782344) × (400 − 295) = 10105.0607335181.
3.2. STEP 2 – calculating volatility. Using the options data selected in the previous step, we approximate
the integrand of the right-hand side of (3.1),(3.2)

z 7→ σ(g(z))2 .

This is attained in three stages: (i) obtaining from the available options a sequence of data points thatwill act as “knots” for the regression based on (3.2), (ii) constructing, using these data points, a functionthat approximates (3.2), and finally (iii) integrating it with respect to the normal density, and obtaining
the implied variance using the identity (3.1). Throughout this step, the bid/ask option prices are used.
STEP 2-1 – converting the data. Each options data consists of the option price and strike price K = Fek ,and this is suﬃcient to compute the corresponding implied volatility σ(k) as defined in (2.4). A data
point (d2 (k), σ(k)2 ) is thereby obtained for each option. More precisely, the forward price F is obtained in
STEP 1, and then the bisection method with error bound 1.0e − 9 is used to compute implied volatilities.
As is shown in Theorem 2.1, d2 (k) is decreasing in k, thereby decreasing in K as well. In order to satisfy
this constraint, we disregard the data points on which this monotonicity fails, and focus on the (longest
possible) interval of data points at which the monotonicity holds. Starting from the put option with thehighest strike (regularly K0 ), each option is visited backward in K to make sure that d2 (k) is decreasing.
Once monotonicity breaks down, the corresponding put option and any put option with lower strikes
are discarded. The same procedure is applied to call options. Starting from the call option with thelowest strike, each option is visited forward in K; once monotonicity breaks down, the associated call
option and any call option with higher strikes are similarly discarded. The first five columns in Table 2
show the option prices selected from Table 1 and their corresponding d2 ’s and implied volatilities. Here,d2 is found to be monotonically decreasing in the strike price, and therefore no truncation is applied.

MODEL-FREE IMPLIED VOLATILITY: FROM SURFACE TO INDEX

Strike
5000
5500
6000
6500
7000
7500
8000
8250
8500
8750
9000
9250
9500
9750
10000
10250
10500
10750
11000
11250
11500
11750
12000
12250
12500
12750
13000
13500
14000
14500

Put
trans. ask bid
1
1
1
1
2
2
1
3
4
3
9
9
4
17
17
16
25
25
20
30
35
30
45
50
45
65
70
65
100 10595
150 150145
215 215205
295 300
295
405 420405
550 570550
730 750730
955
940
1150 1180 1160
1420 1390
1650 1660 1640
1910 1880
2150 2130
2400 2380
2650 2620
2900 2870
3400 3370
3900 3880
4400 4380

Calltrans. ask bid
5120 5090
4620 4590
4120 4090
3620 3590
3120 3100
2630 2600
2150 2140 2110
1890 1870
1650 1630
1420 1390
1190 1170
950 970950
745 765745
565 575
560
400 410400
280 280265
165 175165
100 105
100
55
60
55
30
35
30
18
19
17
9
10
9
6
6
5
3
4
3
2
2
1
1
2
1
1
1
1
1
-

7

diﬀerence type mid price
2133
850
595
350
105
125
385
630
1120
1641
-

Put
Put
Put
Put
Put
Put
Put
Put
Put
Put
Call
Call
Call
Call
Call
Call
Call
Call
Call
-

3.5
16.5
22.5
32.5
47.5
67.5
100
147.5
210
297.5
272.5
170
102.5
57.5
32.5
18
9.5
5.5
3.5
-

Table 1. An example illustrating how the ATM strike and OTM puts and calls are selected.

STEP 2-2 – constructing a function approximating (3.2). After STEP 2-1, we suppose there are M
data points. We denote by x1 ≤ · · · ≤ xM the ascending sequence of d2 (·) and by y1 ≤ · · · ≤ yM thecorresponding sequence of implied volatilities σ(·)2 . In the example in Table 2, M = 19, x1 = −2.3338,

8

M. FUKASAWA, I. ISHIDA, N. MAGHREBI, K. OYA, M. UBUKATA, AND K. YAMAZAKI. . . , xM = 2.322589 and y1 = 0.05886305, . . . , yM = 0.1953966. For each of the following intervals(−∞, x1 ], [x1 , x2 ], [x1 , x2 ], . . . , [xM−1 , xM ] and [xM , ∞),

iv

0.20

0.25

0.30

0.35

0.40

0.45

0.50

a cubic polynomial is defined so that their union is continuous and diﬀerentiable on R and goes througheach data point.
In defining the cubic polynomial for a given interval, say [x j , x j+1 ] for some j = 1, . . . , M − 1, there are
two more “degrees of freedom” remaining after satisfying the requirements that it goes through (x j , y j )and (x j+1 , y j+1 ). It should be noted that there are various ways of using these degrees of freedom. A
typical example uses cubic splines as in Jiang and Tian (2007), which realizes the twice-diﬀerentiability atthese joint points. We choose, however, to obtain polynomials in the way we shall discuss hereafter. For(−∞, x1 ] and [xM , ∞), we use constant extrapolation. In the next section, we will discuss this method incomparison to the one by Jiang and Tian (2007).
{}
The slope of the cubic polynomial at each joint point (x j , y j ); j = 1, . . . , M is first determined. We
initially set y0 (x1 ) = y0 (xM ) = 0. For 2 ≤ j ≤ M − 1, the slope is chosen so that “the incidence angle
becomes equal to the reflected angle”. In other words, we construct a polynomial so that the angle madeby its tangent line at (x j , y j ) and the line going through (x j−1 , y j−1 ) and (x j , y j ) is equal to the angle madeby the same tangent line and the line going through (x j , y j ) and (x j+1 , y j+1 ), as shown in Figure 1. This

−0.20

−0.15

−0.10

−0.05

0.00

0.05

0.10

d2

Figure 1. The incidence angle becomes equal to the reflected angle. The y-axis
refers to IV, which represents Implied Variance in this case.
can be calculated as in the following:)/ ()(
y j+1 − y j y j − y j−1
x j+1 − x j x j − x j−10−−y (x j ) = −
l j+1
lj
l j+1
lj
where l j and l j+1 are Euclidean distances from (x j , y j ) to (x j−1 , y j−1 ) and (x j+1 , y j+1 ), respectively.

MODEL-FREE IMPLIED VOLATILITY: FROM SURFACE TO INDEX

9

For every 1 ≤ j ≤ M, the polynomial on [x j , x j+1 ] is now uniquely determined by the two end-points(x j , y j ) and (x j+1 , y j+1 ) and the slopes at these points. We represent the polynomial byy(x) = a j + b j (x − x j ) + c j (x − x j )2 + d j (x − x j )3 ,x j ≤ x ≤ x j+1

where a j , b j , c j and d j are the coeﬃcients to be determined. It can be easily shown thata j = y j,b j = y0 (x j ),c j = (3∆y j − ∆x j y0 (x j+1 ) − 2∆x j y0 (x j ))/(∆x2j ),d j = (∆y j − y0 (x j )∆x j − c j ∆x2j )/(∆x3j ),where ∆x j = x j+1 − x j and ∆y j = y j+1 − y j . For the polynomials on (−∞, x0 ] and [xM , ∞), we set y = x0
and y = xM , respectively. Figure 2 shows the polynomial based on the example in Table 2.

Figure 2. The cubic polynomial obtained based on the option data in Table 2.

10

M. FUKASAWA, I. ISHIDA, N. MAGHREBI, K. OYA, M. UBUKATA, AND K. YAMAZAKI

strike type priced2a (impl. vol.)
b
c
d
7000
P
3.5 2.322589
0.1953966
0
0
0
8000
P
16.5 1.737578
0.1401579 0.1024657 0.1339089 -0.2523994
8250
P
22.5 1.597871
0.1247173 0.0900612 0.3505619 -1.4609950
8500
P
32.5 1.428667
0.1129279 0.0628586 -0.0399028 0.4739328
8750
P
47.5 1.243389
0.1025435 0.0574971 -0.0524102 0.2406433
9000
P
67.5 1.054255
0.0913947 0.0472180 0.1316943 -0.3684178
9250
P
100.0 0.833485
0.0835569 0.0318685 -0.0201518 0.1658297
9500
P
147.5 0.595460
0.0768361 0.0298054 -0.0284511 0.0918246
9750
P
210.0 0.347682
0.0690620 0.0273430 0.0388834 -0.0912490
10000
P
297.5 0.077152
0.0627555 0.0188023 0.0184341 -0.0065281
10250
C
272.5 -0.211813
0.0586251 0.0146191 -0.0178526 0.0578870
10500
C
170.0 -0.516513
0.0540715 0.0102862 0.0316420 -0.0536746
10750
C
102.5 -0.820640
0.0523597 0.0056111 -0.0151997 0.0501673
11000
C
57.5 -1.128248
0.0506391 0.0020201 0.0231773 -0.0375809
11250
C
32.5 -1.410956
0.0510783 -0.0023874 -0.0067401 0.0342762
11500
C
18.0 -1.678436
0.0519399 -0.0026407 -0.0074597 0.0197729
11750
C
9.5 -1.941339
0.0524815 -0.0067655 0.0380046 -0.0764793
12000
C
5.5 -2.158142
0.0549685 -0.0168207 0.0276429 -0.0136939
12250
C
3.5 -2.333800
0.0588631
0 -0.2828918 0.8919309

Table 2. The values of d2 and a, b, c, d representing the coeﬃcients of the polynomials a j , b j , c j , d j based on the selected option prices listed in Table 1.

3.3. STEP 3 – integrating. Now we obtain the integral of the function. Notice that∫ ∞−∞σ(g(z)) φ(z)dz =2

M−1∑ ∫ x j+1
j=1

σ(g(z)) φ(z)dz +
2

xj

∫ x1
−∞

σ(g(z)) φ(z)dz +2∫ ∞σ(g(z))2 φ(z)dz.

xM

Here, in particular, for every 1 ≤ j ≤ M − 1, we have
∫ x j+1

σ(g(z)) φ(z)dz ≈
2

xj

=

∫ x j+1y(z)φ(z)dz
xj∫ x j+1 ()
a j + b j (z − x j ) + c j (z − x j )2 + d j (z − x j )3 φ(z)dz

xj= a jA j + b jB j + c jC j + d jD j

MODEL-FREE IMPLIED VOLATILITY: FROM SURFACE TO INDEX

11

where
A j := Φ(x j+1 ) − Φ(x j ),()
B j := −(φ(x j+1 ) − φ(x j )) − x j Φ(x j+1 ) − Φ(x j ) ,
(
)
(
)
(
)
C j := − x j+1 φ(x j+1 ) − x j φ(x j ) + 2x j φ(x j+1 ) − φ(x j ) + (1 + x2j ) Φ(x j+1 ) − Φ(x j ) ,D j := (1 − x2j+1 )φ(x j+1 ) − (1 − x2j )φ(x j ) + 3x j (x j+1 φ(x j+1 ) − x j φ(x j ))− 3(1 + x2j )(φ(x j+1 ) − φ(x j )) − x j (3 + x2j )(Φ(x j+1 ) − Φ(x j )).
The implied variance (3.1) can thus be approximated using the above formula.
4. Comparison to the CBOE procedure
This section discusses the new method in comparison to the CBOE procedure. An apparent disadvantage of the CBOE approximation (2.3) is that it leads the underestimation of the expected quadratic
variation due to the reduced contributions of P(K) on K < Kmin and C(K) on K > Kmax . This is exacerbatedeven further when the range of available strike prices is not suﬃciently wide.
Let us illustrate this issue by using the Nikkei 225 options in comparison to the S&P 500 options. On
October 9, 2008, the S&P 500 index dropped to 909.92 USD, and it was followed by the Nikkei stock
average, which decreased to 8276.43 JPY on the next day. With both indexes dropping around 10% on
a single day, the ranges of strikes for the corresponding options are however significantly diﬀerent. As
shown in Table 3, the lowest strike of the Nikkei 225 index option Kmin approaches the spot price of theunderlying index; it is as large as 88% of the spot price, which is tremendously larger than the ratio for
the S&P case. As a result, as can be observed in Table 4, the price of the put option at Kmin becomessignificantly high in comparison to the single-digit prices reported in Table 1 under normal market
conditions. In this case, the cutting-oﬀ of the integral over [0, Kmin ] severely aﬀects the approximation
Spot priceKmin and ratio to spot price on Oct. 9 (Oct. 10)
Oct. 8 (Oct. 9) Oct. 9 (Oct. 10) Nov. maturity
Dec. maturity
S&P
984.94909.92400 (40%)400 (40%)
Nikkei9157.498276.43
7250 (88%)
7250 (88%)

Table 3. The lowest strikes of the S&P 500 and the Nikkei 225 index options andthe spot prices of their underlying market indexes.
because the integrand P(K)/K2 is not suﬃciently small around K = Kmin . Moreover, a sharp decline
in asset value tends to be accompanied with an increase in expected volatility, or equivalently a fatter
tail of P(K). This further amplifies the eﬀect of the cutting-oﬀ. The CBOE procedure (2.3) is therefore
likely to lead to a non-negligible negative bias in options markets with unaccommodative listing rules,
particularly during financial crises. A more accurate estimation of implied volatility is indeed required
for this reason.
Our approach improves this approximation through the process of extrapolation. Jiang and Tian
(2007) employed the cubic spline for the Black-Scholes implied volatility with linear extrapolation.

12

M. FUKASAWA, I. ISHIDA, N. MAGHREBI, K. OYA, M. UBUKATA, AND K. YAMAZAKI

Nov. maturity
Dec. maturity
strike trans. ask bid trans. ask bid
7250
500 520 500 695
695 6557500
580 600 580 605
605 5707750
680 700 680 790
790 7608000
805 805 790 785
900 8558250
925 925 905 930 1020 970
..
..
..
..
..
..
..
.
.
.
.
.
.
.

Table 4. Put option prices of the Nikkei 225 index options with the lowest strikes
on Oct. 10, 2008 maturing in Nov. 2008 and Dec. 2008.

We remark that their method may cause a positive bias due to its inconsistency with the theoretical
asymptotic behavior of the volatility surface described by Lee (2004), where he obtained a model-free√bound σ(k) < 2|k| for large |k|. See also Lemma A.4 for a refinement of this bound. In order to beconsistent with Theorem 2.2 that gives a counterpart of those model-free bounds for our integration
scale, we choose a constant extrapolation as described in STEP 2-2.
Another source of approximation errors in (2.3) is the discretization of the integral with respect to K.
A continuous interpolation would be a more natural treatment since, in theory, P(K) = E[(K − ST )+ ] andC(K) = E[(ST − K)+ ] are continuous functions. They are C1 if ST admits a density. Our algorithm employsa C1 -interpolation by cubic polynomials in the Black-Scholes implied volatility scale. The algorithm by
Jiang and Tian (2007) consists of the C2 -interpolation and numerical integration with the original scale(2.2). However, the C2 -property of these functions is not guaranteed in general, and adhering to this may
cause undesirable oscillations. On the other hand, our C1 -interpolation approach is expected to preventthese oscillations, which is in fact required in the light of Theorem 2.2.
We proceed to integration in the same Black-Scholes implied volatility scale using the formula (1.1).
Since the integral is with the standard normal density, we can utilize the well-known identities for theHermite polynomial system to avoid numerical integration. The model-free formula (1.1) describes a
direct link between the expected quadratic variation and the curvature of the Black-Scholes implied
volatility surface that is familiar to practitioners.
To examine the extent to which our method indeed improves the approximation, let us consider the
Heston model:

(4.1)

√
]
√ [
1
2
2
dSt = St Vt ρdWt + 1 − ρ dWt ,√dVt = λ(v − Vt )dt + η Vt dWt1 ,where (W 1 , W 2 ) is a two-dimensional standard Brownian motion. This model enables us to calculate
easily the theoretical call and put option prices and obtain the following explicit formula for the expected

MODEL-FREE IMPLIED VOLATILITY: FROM SURFACE TO INDEX

13

annualized quadratic variation:
(4.2)

1
1
E[hlog(S)iT ] = E
T
T

[∫ T
0

]
Vt dt = v +1 − e−λT(V0 − v).
λT

See, e.g., Gatheral (2006) for more details. Using artificial options data generated by the Heston pricing
formula, we evaluate the approximations in comparison to the true value of the expected quadratic
variation. We generate call and put option prices with the same set of strikes and nearest two maturities
as the above-mentioned options data on October 10, 2008. The first and second maturities falling in November and December correspond to T = 0.0951864535768645 and T = 0.171898782343988, respectively.
The initial price S0 of the underlying asset on October 10, 2008 is set to 8276.43 JPY and the interest rate
is set to 0. See Table 5 for the generated data for the first maturity. To further approach normal market
conditions, we randomize the bid-ask spread using geometric random variables with success probability
p = 0.8. The ask price is sampled so that it becomes the lowest price above the true value (consistent withthe tick-size used on the OSE) with probability p, and becomes the second lowest price with probability(1 − p)p. It becomes the subsequent prices in the same geometric manner, and the bid price is chosen
similarly. We use the corresponding mid-quote price for each transaction price.
Table 6 shows the expected annualized quadratic variations and the approximated values by the CBOE
procedure and our method based on the following parameter sets:
Parameter Set A: λ = 1, v = 0.2, η = 0.5, ρ = −0.8 and V0 = 0.6;
Parameter Set B: λ = 1, v = 0.2, η = 1.0, ρ = −0.4 and V0 = 0.6;
Parameter Set C: λ = 5, v = 0.04, η = 1.0, ρ = −0.4 and V0 = 0.6;
Parameter Set D: λ = 1.5, v = 0.04, η = 0.3, ρ = −0.7 and V0 = 0.04.
The column “New” stands for the values resulting from the application of our new method. The
parameters in D are typical values obtained by calibrating on a regular day with instantaneous volatility√
V0 = 20%. Under this average volatility level, the tails of the distribution of ST are light, and hence
the approximation errors due to the cutting-oﬀ of tails are negligible. The main source of error lies in
the numerical integration and our method provides better numerical accuracy. The other specifications
A, B and C represent the market conditions at the onset of financial crises, with instantaneous volatility√
V0 ≈ 77%. The results are more striking. The CBOE method severely underestimates the expected
quadratic variation where our method is found to significantly improve the estimation accuracy.

5. New Volatility Index and Empirical evidence
In this section, we evaluate the in-sample forecasting performance of the volatility index from the
Nikkei 225 options based on the new method in comparison to that using the CBOE procedure. The options data were obtained from the Nikkei Economic Electronic Database System (NEEDS) FinancialQuest
database, and contain call and put option prices quoted at the close of trading days from December 22,
1997 to December 1, 2009 (2,935 trading days). Throughout this section, we denote by CSFI-VXJ2 our14

M. FUKASAWA, I. ISHIDA, N. MAGHREBI, K. OYA, M. UBUKATA, AND K. YAMAZAKI

Call
Put
strike ask theoretical bid ask theoretical bid
7250 1370
1361.58 1360 345
335.15 335
7500 1200
1197.41 1190 430
420.98 420
7750 1050
1046.16 1040 525
519.73 515
8000
915
907.99 905 640
631.56 630
8250
790
782.88 780 765
756.45 755
8500
675
670.56 665 910
894.13 890
8750
580
570.59 570 1050
1044.16 1040
9000
485
482.35 475 1210
1205.92 1200
9250
415
405.12 405 1380
1378.69 1370
9500
345
338.08 335 1570
1561.65 1560
9750
290
280.34 280 1760
1753.91 1750
10000 240231.01 230 1960
1954.58 1950
10250 195189.19 185 2170
2162.76 2060
10500 160
153.99 140 2380
2377.56 2370
10750 130124.60 120 2610
2598.17 2590
11000 110100.22 100 2830
2823.79 2820
11250
90
80.15
80 3060
3053.72 3050
11500
70
63.73
60 32903287.30 3280
11750
55
50.40
45 35403523.97 3520
12000
40
39.63
30 37703763.20 3760
12250
40
31.00
30 40104004.57 3900
12500
30
24.12
20 42504247.69 4140
12750
20
18.67
18 4500
4492.24 4290
13000
16
14.38
14 47404737.95 4730
13250
13
11.02
11 49904984.59 4980
13500
10
8.40
8 52405231.97 5230
13750
8
6.38
6 54805479.95 5370
14000
8
4.82
4 57305728.39 5720
14250
5
3.62
3 60005977.19 5970
14500
4
2.71
2 62306226.28 6120
15000
2
1.50
- 67306725.07 6720
15500
1
0.81
- 72407224.38 7220
16000
1
0.44
- 77407724.01 7720
16500
1
0.23
- 82308223.80 8220
17000
1
0.12
- 87308723.69 8720
17500
1
0.06
- 92309223.63 9120

Table 5. Artificial data with maturity November 2008 and parameter set A.

MODEL-FREE IMPLIED VOLATILITY: FROM SURFACE TO INDEXParameter
A
B
C
D

True
0.5816
0.58160.48560.04002008/11
CBOE New
0.4639 0.5767
0.4604 0.5692
0.3930 0.4633
0.0420 0.0408

True
0.5676
0.5676
0.41570.040015

2008/12
CBOE New
0.4041 0.5504
0.4023 0.5460
0.3112 0.4023
0.0378 0.0394

Table 6. Expected annualized quadratic variations.

volatility index and by VXJ2 the index computed based on the CBOE approach.1 Here, the values of
CSFI-VXJ2 and VXJ2 are measured as annualized 30-day quadratic variation, while the volatility indexes
usually refer to their square roots in percentage points.
Table 7 summarizes the descriptive statistics of the CSFI-VXJ2 and VXJ2 volatility indexes. The meanvalue of CSFI-VXJ2 is close to that of VXJ2 , but CSFI-VXJ2 has a slightly larger standard deviation than
VXJ2 . The skewness and kurtosis moments indicate that both CSFI-VXJ2 and VXJ2 are not close to thenormal distribution. The large Ljung-Box statistics up to tenth order (LB(10)) indicate that the null
hypothesis of no autocorrelations is rejected.
For the sake of estimating market variance over the fixed expiration period of 30 days, we employ
ex-post realized variances using high-frequency data. These intraday measures of realized volatility areusually found to provide more accurate variance estimates than those derived from lower-frequency
data such as daily, weekly or monthly observations. The NEEDS historical tick data include the intraday
Nikkei 225 stock average recorded every minute. The Tokyo Stock Exchange (TSE) where the Nikkei
225 index components are traded is open for 9:00-11:00 (morning session) and 12:30-15:00 (afternoon
session). In order to take into consideration the non-trading hours, we employ the adjustment introduced
by Hansen and Lunde (2005), using the past 441 daily returns and realized variances to estimate thetime-varying adjustment parameter.
In order to perform tests that are robust to the microstructure noise, we use as many as 15 estimatorsof realized variance listed in Table 8. A more detailed description of these estimators is provided
in Appendix B. We summarize in Table 9 the descriptive statistics of annualized realized variancescomputed over fixed periods of 30 days, consistent with the expiry of the hypothetical options underlying
the implied volatility indexes. The realized variances have lower means than the CSFI-VXJ2 and VXJ2
indexes, although the standard deviations are rather close to these volatility forecasts. The distribution ofrealized variance is also found to be positively skewed and long-tailed. Judging from LB(10) values, there
are signs of autocorrelations. Figure 3 plots the time-series of market variance along with the CSFI-VXJ2
and VXJ2 indexes. There are sudden surges associated with significant economic and financial shocks
such as the Russian default and Long-Term Capital Management crises in 1998 and the U.S. housing andcredit crisis in 2008.
1

We use these names because the same indexes are computed, updated and made available online to the public
under these titles by the Center for the Study of Finance and Insurance, Osaka University.
http://www-csfi.sigmath.es.osaka-u.ac.jp/structure/activity/vxj.php

16

M. FUKASAWA, I. ISHIDA, N. MAGHREBI, K. OYA, M. UBUKATA, AND K. YAMAZAKI

MeanStandard deviation
Minimum
Maximum
Skewness
Kurtosis
LB(10)

CSFI-VXJ2
0.08160.08240.01200.9454
4.90
31.98
24432

VXJ2
0.08220.07710.01330.8363
4.55
27.31
25271

Table 7. Descriptive statistics of the CSFI-VXJ2 and VXJ2 .
1. Realized variance with returns sampled at the highest frequencies, RV
2. Realized variance with 5-minute returns, RV5min
3. Realized variance with 15-minute returns, RV15min
4. Optimally-sampled realized variance as proposed in Bandi and Russell (2008), RVBR
5. The Bartlett kernel estimator in Barndorﬀ-Nielsen, Hansen, Lunde, and Shephard (2005)
with a finite sample optimal number of autocovariances
BK
proposed by Bandi and Russell (in press), RVBR
6. The two-scale estimator with an asymptotically optimal number of subsamples
proposed by Zhang, Mykland, and Aı̈t-Sahalia (2005), RV ZMA
7. The two-scale estimator with a finite sample optimal number of subsamples
ZMA
proposed by Bandi and Russell (in press), RVBR
8. The bias-corrected two-scale estimator with an asymptotically optimal number ofsubsamples proposed by Zhang et al. (2005), RV BC−ZMA
9. The bias-corrected two-scale estimator with a finite sample optimal number of subsamples
BC−ZMA
proposed by Bandi and Russell (in press), RVBR
10. The flat-top Bartlett kernel estimator with an asymptotically optimal number of
autocovariances proposed by Barndorﬀ-Nielsen, Hansen, Lunde, and Shephard (2008), RV FBK
11. The flat-top Bartlett kernel estimator with a finite sample optimal number of
FBK
autocovariances proposed by Bandi and Russell (in press), RVBR
12. The flat-top cubic kernel estimator with an asymptotically optimal number of
autocovariances proposed by Barndorﬀ-Nielsen et al. (2008), RV FCK
13. The flat-top cubic kernel estimator with a finite sample optimal number ofFCK
autocovariances proposed by Bandi and Russell (in press), RVBR
14. The flat-top modified Tukey-Hanning kernel estimator with an asymptotically
optimal number of autocovariances proposed by Barndorﬀ-Nielsen et al. (2008), RV FMTH
15. The flat-top modified Tukey-Hanning kernel estimator with a finite sample optimal
FMTH
number of autocovariances proposed by Bandi and Russell (in press), RVBR

Table 8. List of estimators of realized variance using high-frequency data.

MODEL-FREE IMPLIED VOLATILITY: FROM SURFACE TO INDEXRV
RV5min
RV15min
RVBR
BK
RVBR
RV ZMAZMARVBR
RV BC−ZMA
BC−ZMARVBR
RV FBKFBKRVBR
RV FCKFCKRVBR
RV FMTHFMTH
RVBR

Mean Standard error Minimum Maximum Skewness Kurtosis LB(10)0.06400.08060.00780.9334
7.10
61.12
25977
0.06280.08190.00630.9346
7.03
60.09
25914
0.06240.09280.00671.0650
7.21
61.30
26031
0.06320.08700.00611.0066
7.37
65.08
25750
0.06280.08960.0064
1.0577
7.55
67.60
25685
0.06330.08710.00671.0265
7.30
63.67
25889
0.06300.09070.0061
1.0684
7.50
66.92
25689
0.06430.08070.00640.9305
7.06
61.48
25695
0.06340.08510.00620.9865
7.28
64.26
25663
0.06360.08340.00620.9601
7.15
62.43
25690
0.06360.08450.0062
0.9762
7.23
63.51
25672
0.06360.08600.0061
0.9946
7.28
64.15
25664
0.06350.08410.0062
0.9716
7.21
63.33
25670
0.06360.08560.00620.9888
7.23
63.43
25697
0.06360.08410.00620.9694
7.15
62.41
25700

Table 9. Descriptive statistics of realized variances.
1

Market variance
NVXJ
VXJ

0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

20090518

20080717

20070921

2006113020060210

20050418

20040624

20030828

20021106

20020117

20010326

20000602

19990809

19981014

19971222

Figure 3. Time-series of realized variance and implied variance indexes.

17

18

M. FUKASAWA, I. ISHIDA, N. MAGHREBI, K. OYA, M. UBUKATA, AND K. YAMAZAKI

A formal evaluation of forecasting performance of the various implied volatility indexes can be madeusing the following regression modelV̂t = α + β IVt + ut ,(5.1)

where V̂ and IV represent the annualized realized variance and implied variance, respectively. Table 10
reports the coeﬃcient of determination R2 from the regression equation (5.1). Irrespective of the estimatorof realized variance used, the R2 values obtained using CSFI-VXJ2 are found to be higher than thoseobtained using VXJ2 . These results show that the CSFI-VXJ2 forecasts of return variability are moreaccurate than those based on VXJ2 .
Recall that the indexes are constructed as approximations of the expected quadratic variation. Although the expectation is with respect to a risk-neutral measure, by assuming that it coincides withthe physical measure, we further assess the forecast error of CSFI-VXJ2 and VXJ2 using four loss functions: the root mean squared error (RMSE), the root mean squared percentage error (RMSPE), the mean
absolute error (MAE) and the mean absolute percentage error (MAPE) defined, respectively, as
v
u
t
N
√
∑
1 ∑ ( IV − V̂ )2
N
1
2
RMSE = N i=1 (IV − V̂) ,,RMSPE =
N
V̂
i=1

MAE

= N1

∑N

i=1 | IV − V̂ |,

MAPE

=

1
N

N
∑
i=1

IV − V̂V̂.

Here we use N = 2, 935 for the number of observations of realized variances and implied variance indexes
over the sample period. Table 11 displays the estimates of biases and the loss functions for CSFI-VXJ2
and VXJ2 . The bias is always positive, but it is smaller in magnitude for CSFI-VXJ2 than for VXJ2 . The
values of loss functions for CSFI-VXJ2 are uniformly smaller than those for VXJ2 . Judging from theseresults, the volatility index proposed in this paper provides a better approximation of short-term market
volatility and has a better forecasting performance than the index computed using the CBOE procedure.
6. Conclusion
This study proposes a new model-free approach to approximating the expected quadratic variations ofasset prices based on related options premia. The relation between quadratic variations and option prices
is model-free and it is illustrated within the framework of the familiar Black-Scholes implied volatility
scale. The new approximation method avoids numerical integration by taking advantage of the formof the integral with respect to the standard normal density and using the polynomial interpolation of
the integrand. Based on this approach, a new volatility index is developed using the Nikkei 225 index
options. The new volatility benchmark provides good approximations of the true values generated underthe Heston stochastic volatility model. This is indicative of better numerical eﬃciency than alternative
approaches such as the CBOE procedure, which rely on discretization to evaluate integrals. Basedon various estimators of realized variance, the empirical evidence from market prices also suggeststhat the new volatility index is associated with better forecast accuracy. The model-free formula forthe expected quadratic variations is thus conducive to the reduction of approximation errors and to
improved numerical eﬃciency and forecasting accuracy.

MODEL-FREE IMPLIED VOLATILITY: FROM SURFACE TO INDEX

market variance
RV
RV5min
RV15min
RVBR
BK
RVBR
RV ZMAZMARVBR
RV BC−ZMA
BC−ZMARVBR
RV FBKFBKRVBR
RV FCKFCKRVBR
RV FMTHFMTH
RVBR

CSFI-VXJ2
0.488
0.501
0.427
0.482
0.460
0.448
0.458
0.504
0.491
0.501
0.496
0.490
0.497
0.490
0.498

19

VXJ2
0.457
0.471
0.398
0.451
0.430
0.417
0.427
0.473
0.459
0.469
0.464
0.458
0.465
0.459
0.466

Table 10. Coeﬃcients of determination R2 from regressions of realized varianceon implied variance indexes.
Acknowledgements
The authors are thankful for the financial support from the Osaka Securities Exchange and the Japan
Science and Technology Agency CREST. Masaaki Fukasawa and Kazutoshi Yamazaki acknowledge thegrants-in-aid for scientific research No. 21740074 and No. 22710143, respectively, from the Japanese
Ministry for Education, Culture, Sports, Science and Technology. Kosuke Oya and Masato Ubukata arethankful for the grants-in-aid for scientific research No. 22243021 from the Japanese Society for the Promotion of Science. Isao Ishida also acknowledges the grants-in-aid for scientific research No. 20530265.
Nabil Maghrebi is thankful for the grants-in-aid No. 21330077, and the Center for the Study of Financeand Insurance where this study was completed while he was Visiting Professor. The authors also wish
to thank the Editor and anonymous referee for useful comments and suggestions from which this paper
has greatly benefited.
Appendix A. Pricing formula for log contract
In this appendix, we give the proofs of Theorems 2.1 and 2.2. At first, let us observe that theBlack-Scholes implied volatility σ : R → [0, ∞) is well-defined. The value PBS (k, σ) coincides with theundiscounted version of the Black-Scholes put option price with strike K = Fek , maturity T and volatility
σ > 0. It is well-known that the Black-Scholes price is an increasing function of the volatility parameter.
Hence the inverse function PBS (k, ·)−1 is well-defined. Besides, it holds thatlim PBS (k, σ) = F(ek − 1)+

σ→0

and

lim PBS (k, σ) = Fek .

σ→∞Table 11. Estimates of bias and four loss functions for the CSFI-VXJ2 and VXJ2 indexes.

market
bias
RMSE
RMSER
MAE
MAPE
2
2
2
2
2
2
2
2
variance
CSFI-VXJVXJCSFI-VXJVXJCSFI-VXJVXJCSFI-VXJVXJCSFI-VXJ2 VXJ2
RV
0.01750.01820.0657
0.06610.68320.73950.03190.03310.51710.5673RV5min0.01880.01940.06550.06600.73830.80150.03220.03350.55480.6064
RV15min0.01920.01980.07600.07690.84740.9086
0.03630.03770.65100.7051
RVBR
0.01840.01900.06880.06970.75410.81920.03290.03430.5656
0.6188BK
RVBR
0.01870.01930.07180.07270.74260.80040.0336
0.03490.56690.6172RV ZMA
0.01830.01890.07150.07220.75970.82700.03400.03520.57710.6309ZMA
RVBR
0.01850.01920.07240.07340.76150.8226
0.03400.03530.5812
0.6328BC−ZMA
RV
0.01730.01790.06450.06490.70920.77850.03120.03240.52410.5773BC−ZMA
RVBR
0.01810.01880.06740.06820.72730.79290.03240.03370.54900.6012
FBK
RV
0.01800.0186
0.06590.06660.72490.79250.03200.03330.54590.5985FBK
RVBR
0.01800.0186
0.06670.06740.72730.79440.03220.03350.54800.6003
FCK
RV
0.01800.0186
0.06770.06860.73350.79960.03250.03380.55360.6059FCK
RVBR
0.01800.0186
0.06650.06720.72600.79330.03210.03340.54730.5996FMTH
RV
0.01800.0186
0.06750.06840.73720.80320.0326
0.03390.55660.6090FMTH
RVBR
0.01800.0186
0.06640.06710.73090.79870.03230.0336
0.55210.6047

20
M. FUKASAWA, I. ISHIDA, N. MAGHREBI, K. OYA, M. UBUKATA, AND K. YAMAZAKI

MODEL-FREE IMPLIED VOLATILITY: FROM SURFACE TO INDEX

21

Because(K − F)+ ≤ P(K) < K
for all K > 0 due to Jensen’s inequality, σ(k) is well-defined for all k ∈ R.
Since ST admits a density by assumption, we have
D(K) := E[K > ST ] =

dP
(K).
dK

Notice that D is an increasing function of K. Let us suppose, for simplicity, that T = 1 hereafter. Let
DBS (K) be a function of K = Fek defined as
DBS (K) =dPBS(log(K/F), Σ)|Σ=σ(log(K/F)) = Φ(−d2 (k, σ(k))).
dK

By definition,
dP
(K)
dK
d
=
PBS (log(K/F), σ(log(K/F)))dK
1 ∂PBSdσ= DBS (K) +(log(K/F), σ(log(K/F))) (log(K/F))
K ∂Σ
dk
dσ
= DBS (K) + φ(−d2 (log(K/F), σ(log(K/F)))) (log(K/F)).
dk

D(K) =(A.1)

Lemma A.1. For every k ∈ R, we have−d2 (k, σ(k))dσ(k) < 1.
dk

Proof: Put f (k) := −d2 (k, σ(k)). The inequality is trivial when f (k) = 0. If f (k) > 0, it follows from (A.1)
that
1 − Φ( f (k))
D(Fek ) − DBS (Fek )dσf (k) (k) = f (k)≤ f (k)< 1.
dk
φ( f (k))φ( f (k))
Here we used the fact that 0 ≤ D(K) ≤ 1 by definition and a well-known estimate
1 − Φ(x) < x−1 φ(x), x > 0.

(A.2)
For the case f (k) < 0, we havef (k)Φ( f (k))
1 − Φ(− f (k))
D(Fek ) − DBS (Fek )dσ(k) = f (k)≤ − f (k)= − f (k)< 1.
dk
φ( f (k))φ( f (k))φ(− f (k))////

Lemma A.2. For all k with d1 (k, σ(k)) ≥ 0, we have−d1 (k, σ(k))dσ(k) < 1.
dk

22

M. FUKASAWA, I. ISHIDA, N. MAGHREBI, K. OYA, M. UBUKATA, AND K. YAMAZAKI

Proof: By definition, it holds that for all K > 0,
KD(K) ≥ P(K).
Combining this and (A.1), we haveFΦ(−d1 (k, σ(k))) + Kφ(−d2 (k, σ(k)))dσ(k) ≥ 0dk

with k = log(K/F). Since Kφ(−d2 ) = Fφ(−d1 ), we obtain by (A.2)
1 − Φ(d1 (k, σ(k)))1dσ(k) ≥ −>−.
dk
φ(d1 (k, σ(k)))
d1 (k, σ(k))////
Lemma A.3. The mapping k 7→ −d2 (k, σ(k)) is increasing.
Proof: Put f (k) := −d2 (k, σ(k)). By definition,f (k) =and therefore
(A.3)

σ(k)
k
+
,
σ(k)
2

{
}
{
}
df
1
dσ
k
1 dσ
1
dσ
dσ
(k) =
1−
(k)
+
(k) =
1−
(k) f (k) +
(k).
dk
σ(k)
dk σ(k)
2 dk
σ(k)
dk
dk

Hence, by Lemma A.1, we have
df
dσ
(k) >
(k).
dk
dk
It suﬃces then to treat the case dσ/dk < 0. By rewriting (A.3), we have
{
}
df
1
dσ
(k) =
1+
(k)d1 (k, σ(k)) .
dk
σ(k)
dk
If d1 (k, σ(k)) < 0, we have d f /dk > 0 under dσ/dk < 0. If d1 (k, σ(k)) ≥ 0, we can use Lemma A.2 to obtainthe same inequality.
////
Lemma A.4. It holds for k ≥ 0 that−d2 (k, σ(k)) ≥

√
2k

and

dσ
1
(k) < √ .
dk
2k√
2k∗ and it holds for all k > k0 ≥ k∗ that√√√
σ(k) < σ(k0 ) − 2k0 + 2k ≤ 2k.

Moreover, there exists k∗ > 0 such that σ(k∗ ) =Proof: The first inequality comes from the fact that the arithmetic mean exceeds the geometric mean.
The second follows from this and Lemma A.1. To consider the last claim, observe that
E[(ST − K)+ ] = E[(K − ST )+ ] + F − K = FΦ(d1 (k, σ(k))) − KΦ(d2 (k, σ(k)))
by definition, or “Call-Put Parity”. The left-hand side goes to 0 as K → ∞ and by (A.2)
KΦ(d2 (k, σ(k))) = Fek φ(−d2 (k, σ(k)))

1 − Φ(−d2 (k, σ(k)))
F< √ → 0 as k → ∞.
φ(−d2 (k, σ(k))))2 πk

MODEL-FREE IMPLIED VOLATILITY: FROM SURFACE TO INDEX

23√
Hence d1 (k) → −∞ as k → ∞, which implies that σ(k) < 2k for suﬃciently large k. The last inequality√
holds since 2k − σ(k) is increasing by the second inequality.
////
Here, notice that results similar to the previous lemma were obtained by Lee (2004) and Rogers andTehranchi (2010).
−p

Lemma A.5. If there exists p > 0 such that ST is integrable, then we have(A.4)

lim σ(k)φ(d2 (k, σ(k))) = 0

and

k→±∞

lim k

k→±∞

dσ
(k)φ(d2 (k, σ(k))) = 0.
dk

Proof: By the results of Lee (2004), there exist k∗ < 0 and β ∈ (0, 2) such that
−

σ(k)2
<β
k

for all k < k∗ . It follows thatσ(k)kd2 (k, σ(k)) = −
−
>
σ(k)
2

√

|k|
−
β

√
β|k| 2 − β √
= √
|k|.
2
2 β

With the aid of Lemma A.4, these estimates imply that φ(d2 (k, σ(k))) decays exponentially fast as |k| → ∞and that the first claim of (A.4) holds. To see the second claim, notice that by (A.1)
D(Fek ) − DBS (Fek ) = φ(−d2 (k, σ(k)))

dσ
(k),
dk

so that it now suﬃces to show
lim kD(Fek ) = 0, lim k(1 − D(Fek )) = 0,

k→−∞

k→∞

lim kDBS (Fe ) = 0, lim k(1 − DBS (Fek )) = 0.
k

k→−∞

k→∞

The first-line equations follow from the integrability condition of ST , whereas the remaining two equations are obvious from the fact that DBS (Fek ) = Φ(−d2 (k, σ(k))).
////Proof of Theorem 2.1: Since the second derivative of P coincides with the density of ST ,∫ ∞d2 P
−2E[log(ST /F)] = −2log(K/F) 2 (K)dK
dK
0
∫ ∞ 2
d P
= −2
k 2 (Fek )Fek dk.
−∞ dK
By Lemma A.3, the inverse function g of d2 is well-defined. Using (A.1), we haved2 P kdD k(Fe )(Fe ) =
2
dK
dK
{
}
)
dg−1 (
1
dσ
d2 σ
−1
−1
= k φ(g (k)) −(k) 1 + g (k) (k) + 2 (k)
dk
dk
dk
Fe
Since

dg−1
d
φ(g−1 (k)) = −φ(g−1 (k))g−1 (k)
(k),
dk
dk

24

M. FUKASAWA, I. ISHIDA, N. MAGHREBI, K. OYA, M. UBUKATA, AND K. YAMAZAKI

we have−∫ ∞By (A.4), we obtain

Since

∫ ∞

−∞

kφ(g−1 (k))g−1 (k)]∞[dg−1 dσdσ(k) (k)dk = k (k)φ(g−1 (k))
dk
dk
dk
−∞
}
∫ ∞{
dσ
d2 σ
−
(k) + k 2 (k) φ(g−1 (k))dk.
dk
dk
−∞

{
}
dg−1
dσ
−2E[log(ST /F)] = 2
φ(g (k)) k(k) +(k) dk.
dk
dk
−∞
∫ ∞

−1

∫ ∞
[
]∞
dg−1
dσ
−1
(k)dk
φ(g (k)) (k)dk = φ(g (k))σ(k)+φ(g−1 (k))g−1 (k)σ(k)
−∞
dk
dk
−∞
−∞
−1

and
k + g−1 (k)σ(k) = k + d2 (k, σ(k))σ(k) = −by definition, we obtain−2E[log(ST /F)] = −

∫ ∞
−∞

−1

φ(g (k))σ(k)

−1
2 dg

dk

(k)dk =

σ(k)2
2
∫ ∞
−∞

σ(g(z))2 φ(z)dz.

Here we have used (A.4).
////
Proof of Theorem 2.2: The first inequality for the derivative of σ̂ follows from (A.3) and Lemma A.1.
The existence of z∗ is equivalent to the existence of k∗ in Lemma A.4. The other inequalities follow fromthe fact thatk 7→ d2 (k, σ(k))2 − d1 (k, σ(k))2 = 2k////

is an increasing function.
Appendix B. Realized variance estimators

Here, we start with a brief review of various estimators of realized variance employed in Section
5. It is assumed that the logarithmic equilibrium price follows a continuous semimartingale processdp∗ (s) = µ(s)ds + σ(s)dW(s), where µ(s) is a drift coeﬃcient, W(s) is a standard Brownian motion, and σ2 (s)is the instantaneous variance of the equilibrium price p∗ . We consider that the parameter of interest is∫ t+1
an integrated variance over the t-th trading day, [t, t + 1), defined as Vt := t σ2t dt. Throughout thissection, let E be the expectation operator under the physical measure.
B.1. RV: Realized variance with returns sampled at the highest frequencies. Let the i-th intraday priceon a day t be denoted by p(ti ). One way to estimate the integrated variance V is to use the sum of squared
returns:
n (
∑
)2
(B.1)
p(ti ) − p(ti−1 ) ,
RV =
i=1

where n represents the number of observed intraday returns in [t, t + 1). The estimator is generally called
realized variance or realized volatility. If p(ti ) is equal to the equilibrium price p∗ (ti ), then RV provides aconsistent estimate of V.

MODEL-FREE IMPLIED VOLATILITY: FROM SURFACE TO INDEX

25

However, RV fails to satisfy the consistency condition when there is market microstructure noise as
usually documented in real high-frequency data. The microstructure noise can be induced by various
market frictions such as the discreteness of price changes, bid-ask bounces, and asymmetric information
across traders, inter alia. For the rest of this appendix, assume the i-th intraday return ri is contaminatedby microstructure noise as follows:
ri = p∗ (ti ) − p∗ (ti−1 ) + η(ti ) − η(ti−1 )|{z} | {z }

(B.2)

r∗i

i

where η represents microstructure noise. A growing literature attempts to examine diﬀerent estimatorsof integrated variance from microstructure noise-contaminated high-frequency data, including Zhang
et.al. (2005), inter alia.
B.2. RV5min and RV15min : Realized variance with 5- and 15-minute returns. Another classical method
uses the realized variance constructed from intraday returns sampled at moderate frequencies rather
than at the highest frequency because the realized variance with finer frequency is more sensitive to
microstructure noise. This approach can partially oﬀset the bias induced by the microstructure eﬀects. In
practice, researchers are necessarily forced to select a moderate sampling frequency. For example, it may
be regarded as around those frequencies for which realized variance signature plots under alternative
sampling frequencies are leveled oﬀ. Evidence from previous studies suggests that it is optimal to use 5
to 30-minute return data. Hence, we employ RV5min and RV15min which are equal to the sum of squared
5- and 15-minute returns.
B.3. RVBR : Optimally-sampled realized variance. There is a trade-oﬀ between the microstructure noiseinduced bias and variance reduction at high sampling frequencies. To take this trade-oﬀ into account,
Bandi and Russell (2008) provide a theoretical justification for the choice of optimal sampling frequency
based on the mean squared error (MSE) criterion. They derive the following approximation of the
) 13
(
IQ
∗
∗
optimal number of observations n based on MSE minimization in a finite sample n ≈ (E(2 ))2 , where()∫T
IQ represents an integrated quarticity of the equilibrium price process IQ = 0 σ4 (s)ds . It is estimated∑ˆ = n n r4 with low frequency returns such as 15-minute returns. Following
by realized quarticity IQ
3

i=1 i

the consistent estimator of noise moment as shown by Bandi and Russell (2008), E(2 ) can be estimated(()2 )1/3∑
ˆ Ê(2 )by Ê(2 ) = 1 n r2 at the highest frequency. Thus, RVBR is equal to RV with n̂∗ = IQ/
.
n

i=1 i

BK : The Bartlett-type kernel estimator in Barndorﬀ-Nielsen et al. (2005) with a finite sample
B.4. RVBR
optimal number of autocovariances proposed by Bandi and Russell (in press). RV5min , RV15min and
RVBR have the obvious drawback that they do not incorporate all available observations, and useful
information may thereby be lost. The problem of estimating the integrated variance under microstructurenoise is similar to the autocorrelation corrections that are used in the estimation of long-run variance instationary time-series (see, for instance, Newey and West (1987) and Andrews (1991)). So it is natural toconsider kernel-based estimators of integrated variance under microstructure noise. Barndorﬀ-Nielsen

26

M. FUKASAWA, I. ISHIDA, N. MAGHREBI, K. OYA, M. UBUKATA, AND K. YAMAZAKI

et al. (2005) examine the Bartlett-type kernel estimator defined as(B.3)

RV BK =

H (∑(n − 1 H − 1)
H − h)
γ0 + 2
γh ,
n
H
H
h=1

∑n−h
where γh = i=1
ri ri+h is the h-th autocovariance of intraday returns and γ0 is equal to realized varianceusing returns sampled at the highest frequencies. This estimator weights the realized variance and theH-th return autocovariances by Bartlett weights. The optimal number of autocovariances is obtained
through MSE minimization in finite sample (see equations 7 to 10 in Bandi and Russell (in press) for theexact MSE minimization expressions). There is a convenient rule-of-thumb for choosing H in practice) 31(
as proposed in Bandi and Russell (in press). The expression is obtained as HBR ≈ 2n3V
n. V and2 IQ
IQ are estimated using realized variance and realized quarticity with lower frequency returns such as
15-minute returns. Hence, RV BK with finite-sample optimal number of autocovariances HBR leads to
BK .
RVBR
B.5. RV ZMA : The two-scale estimator with an asymptotically optimal number of subsamples proposedby Zhang et al. (2005). A two-scale or subsampling estimator is proposed by Zhang et al. (2005) in thespirit of the estimation of the long-run variance studied by Carlstein (1986). Denote the original grid ofobservation times as Ψ = {t0 , t1 , t2 . . . , tn }. Consider Ψ is partitioned into K̃ nonoverlapping subgrids,(j)(j)
Ψ = {t j−1 , t j−1+K̃ , t j−1+2K̃ , . . . }, j = 1, . . . , K̃, Then, the realized variance for the subgrid Ψ is defined as
K̃

(B.4)

K̃

(j)
RV =
K̃

nj
∑
(
)2
p(t(j−1)+iK̃ ) − p(t(j−1)+(i−1)K̃ ) .
i=1

The two-scale estimator in Zhang et al. (2005) is given by(B.5)

ZMA = (1/K̃)

K̃
∑
j=1

RV

( j)
− (n̄/n)RV,
K̃

where n̄ = (n − K̃ + 1)/K̃ and RV is the realized variance for the full grid Ψ. The second term corrects
the bias in the first term. The asymptotic optimal number of subsamples K̃ZMA derived by minimizing( 3(E(2 ))2 )1/3
the estimator’s asymptotic variance is given by K̃ZMA =n2/3 . IQ and E(2 ) are estimated by
IQ∑
realized quarticity with 15-minute returns and Ê(2 ) = n1 ni=1 r2i at the highest frequency, respectively.
Thus, ZMA with K̃ZMA leads to RV ZMA .
ZMA
B.6. RVBR: The two-scale estimator with a finite-sample optimal number of subsamples proposedby Bandi and Russell (in press). Barndorﬀ-Nielsen et al. (2005) show that the two-scale estimator isalmost identical to the modified Bartlett kernel estimator. Bandi and Russell (in press) additionally showthat the finite sample MSEs of RV BK and ZMA are very similar in practice. Hence, the ZMA with K̃ = HBRZMAis represented by RVBR.

MODEL-FREE IMPLIED VOLATILITY: FROM SURFACE TO INDEX

27

B.7. RV BC−ZMA : The bias-corrected two-scale estimator with an asymptotically optimal number ofsubsamples proposed by Zhang et al. (2005). The two-scale estimator ZMA has a finite-sample bias asshown in Zhang et al. (2005) who provide an approximate correction for this bias. On the other hand,Bandi and Russell (in press) report the exact bias-correction form. The bias-corrected estimator can be
defined, following a suggestion by Bandi and Russell (in press), as( K̃n − 1 + 2K̃ − K̃2 − n )−1(B.6)
BC(ZMA) = c(K̃, n)ZMA, c(K̃, n) =.
K̃nSince BC(ZMA) is asymptotically equivalent to ZMA, the asymptotically optimal number of subsamples
is given by K̃ZMA . Thus, BC(ZMA) with K̃ZMA can be described by RV BC−ZMA .
BC−ZMA: The bias-corrected two-scale estimator with a finite-sample optimal number of subB.8. RVBR
samples proposed by Bandi and Russell (in press). Since BC(ZMA) is unbiased in a finite-sample, theoptimal number of subsamples is provided by minimizing the finite-sample variance of BC(ZMA). Bandiand Russell (2008, in press) show that the optimal number of subsamples is defined as
[
[(
]
(
)]
)2
K̃BR = arg min Var BC(ZMA) = arg min c(K, n) Var(ZMA) ,(B.7)
0<K̃/n≤1/2

0<K̃/n≤1/2

where Var(ZMA) can be expressed in equation 6 of Bandi and Russell (2008). Hence, BC(ZMA) with K̃BR
BC−ZMAcan be represented by RVBR.
B.9. RV FBK : The flat-top Bartlett kernel estimator with an asymptotically optimal number of autocovariances proposed by Barndorﬀ-Nielsen, Hansen, Lunde, and Shephard (2008). Barndorﬀ-Nielsen etal. (2008) examine the following unbiased flat-top kernel type estimator (called realized kernel)(B.8)

RK = γ0 +H∑k(x)(γh + γ−h ),h=1∑nwhere the non-stochastic k(x) ∈ [0, 1] for x = h−1i=1 ri ri−h with h =H is a weight function and γh =−H, · · · , H. The flat-top Bartlett kernel estimator is equivalent to RK in case where k(x) = 1 − x. For
this class of kernels, Barndorﬀ-Nielsen et al. (2008) show that the asymptotic distribution of RK − V is
mixed normal with zero mean and rate of convergence n1/6 when H = cn2/3 where c is a constant. Then,4
the asymptotically optimal value of c which minimizes the asymptotic variance is given by c∗ ≈ 2.28ξ 3 ,√where ξ2 = σ2η / IQ. Hence, RK with k(x) = 1 − x and H = c∗ n2/3 corresponds to RV FBK .
B.10. RV FCK and RV FMTH : The flat-top cubic kernel estimator and the flat-top modified Tukey-Hanningkernel estimator with an asymptotically optimal number of autocovariances proposed by BarndorﬀNielsen et al. (2008). The estimators based on the cubic kernel and the modified Tukey-Hanning kernel
are equivalent to RK with k(x) = 1 − 3x2 + 2x3 and k(x) = {1 − cosπ(1 − x)2 }/2. When H = cn1/2 , RK for thisclass of kernels is consistent at the rate of convergence n1/4 as shown in Barndorﬀ-Nielsen et al. (2008).
The asymptotically optimal value of c is expressed as
v
u
√
t
}
1,1 {
k
3k•0.0 k•2,2•∗
c = ρ 0,0 1 +
(B.9)
,
k•
ρ(k•1,1 )2

28

M. FUKASAWA, I. ISHIDA, N. MAGHREBI, K. OYA, M. UBUKATA, AND K. YAMAZAKI

∫1
∫1
∫1
√
where ρ = V/ IQ, k•0,0 = 0 k(x)2 dx, k•1,1 = 0 k0 (x)2 dx and k•2,2 = 0 k00 (x)2 dx, where the primes represent

derivatives. The values of (k•0,0 , k•1,1 , k•2,2 ) amount to (k•0,0 , k•1,1 , k•2,2 ) = (0.371, 1.20, 12.0) for cubic kernel and(k•0,0 , k•1,1 , k•2,2 ) = (0.219, 1.71, 41.7) for modified Tukey-Hanning kernel. We define RV FCK and RV FMTH as
RK with H = c∗ n1/2 at k(x) = 1 − 3x2 + 2x3 and k(x) = {1 − cosπ(1 − x)2 }/2.

FBK , RV FCK and RV FMTH : The flat-top Bartlett kernel estimator, the flat-top cubic kernel
B.11. RVBR
BR
BR
estimator and the flat-top modified Tukey-Hanning kernel estimator with a finite-sample optimal
number of autocovariances proposed by Bandi and Russell (in press). Bandi and Russell (in press)
provide an alternative way to choose the number of autocovariances in finite samples. Denote H as φnwith 0 < φ ≤ 1. The optimal value of φ is defined as[](B.10)
φ∗ = arg min (bias(RK))2 + Var(RK) ,0<φ≤1

where bias(RK) = 0 and Var(RK) can be expressed in Theorem 3 of Bandi and Russell (in press). Thus,FBK ,RK with H = φ∗ n for Bartlett kernel, cubic kernel and modified Tukey-Hanning kernel leads to RVBR
FCK
FMTH
RVBR
and RVBR, respectively.
References[1] Andrews, D. W. K. (1991): Heteroskedasticity and autocorrelation consistent covariance matrix estimation.
Econometrica 59, 817-858.
[2] Bandi, F. M. and Russell, J.R. (2008): Microstructure noise, realized variance, and optimal sampling. Review of
Economic Studies 75, 339-369.
[3] Bandi, F. M. and Russell, J.R. (in press): Market microstructure noise, integrated variance estimators, and theaccuracy of asymptotic approximations. Journal of Econometrics.
[4] Barndorﬀ-Nielsen, O.E., Hansen, P.R., Lunde, A., and Shephard N. (2008): Designing realized kernels tomeasure the ex-post variation of equity prices in the presence of noise. Econometrica 76, 1481-1536.
[5] Carlstein, E. (1986): The use of subseries values for estimating the variance of a general statistic from astationary sequence. Annals of Statistics 14, 1171-1179.
[6] Carr, P. and Lee, R. (2009): Volatility Derivatives, Annu. Rev. Financ. Econ. 1:1-21.
[7] Carr, P., Lee, R. and Wu, L. (2009): Variance swaps on time-changed Lévy processes. Finance and Stochastics,
forthcoming.
[8] Carr, P. and Wu, L. (2006): A tale of two indices. Journal of Derivatives 13(3), 13-29.
[9] Chriss, N. and Morokoﬀ, W. (1999): Market Risk for Volatility and Variance swaps, Risk, October 1999.
[10] Fukasawa, M. (2010): The normalizing transformation of the implied volatility smile, Mathematical Finance,
forthcoming.
[11] Gatheral, J. (2006): The volatility surface; a practitioner’s guide, Wiley.
[12] Hansen, P. R. and Lunde, A. (2005): A realized variance for the whole day based on intermittent highfrequency data. Journal of Financial Econometrics 3, 525-554.
[13] Heston, S.L. (1993): A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bondand Currency Options. The Review of Financial Studies Vol. 6, No. 2, pp. 327-343.
[14] Jiang, G.T and Tian, Y.S. (2007): Extracting Model-Free Volatility from Option Prices: An Examination of the
VIX Index, Journal of Derivatives, 35-60

MODEL-FREE IMPLIED VOLATILITY: FROM SURFACE TO INDEX

29[15] Lee, R. W. (2004): The moment formula for implied volatility at extreme strikes. Math. Finance 14, no. 3,
469–480.
[16] Morokoﬀ, W., Akesson, Y. and Zhou, Y. (1999): Risk management of volatility and variance swaps, Firmwide
Risk Quantitative Modeling Notes, Goldman Sachs & Co.
[17] Newey, W.K. and West, K.D. (1987): A simple, positive semi-definite, heteroskedasticity and autocorrelationconsistent covariance matrix. Econometrica 55, 703-708.
[18] Rogers, C. and Tehranchi, M. (2010): Can the implied volatility surface move by parallel shifts? Finance andStochastics 14(2), 235–248. (2010)[19] Zhang, L., Mykland, P.A., and Aı̈t-Sahalia, Y. (2005): A tale of two time scales: determining integrated
volatility with noisy high-frequency data. Journal of the American Statistical Association 100, 1394-1411.

(Masaaki Fukasawa) Center for the Study of Finance and Insurance, Osaka University, 1-3 Machikaneyama,
Toyonaka, Osaka, 560-8531 Japan.
E-mail address: fukasawa@sigmath.es.osaka-u.ac.jp(Isao Ishida) Center for the Study of Finance and Insurance, Osaka University, 1-3 Machikaneyama, Toyonaka, Osaka, 560-8531 Japan.
E-mail address: i-ishida@sigmath.es.osaka-u.ac.jp(Nabil Maghrebi) Graduate School of Economics, Wakayama University, 930 Sakaedani, Wakayama, 640-8510
Japan
E-mail address: nebilmg@eco.wakayama-u.ac.jp(Kosuke Oya) Graduate School of Economics, Osaka University, 1-7, Machikaneyama Toyonaka, Osaka,
560-0043 Japan
E-mail address: oya@econ.osaka-u.ac.jp(Masato Ubukata) Department of Economics, Kushiro Public University of Economics, 4-1-1 Ashino, Kushiro,
Hokkaido, 085-8585 Japan
E-mail address: ubukata@kushiro-pu.ac.jp(Kazutoshi Yamazaki) Center for the Study of Finance and Insurance, Osaka University,1-3 Machikaneyama,
Toyonaka, Osaka, 560-8531 Japan.
E-mail address: k-yamazaki@sigmath.es.osaka-u.ac.jp