text:
\section{Review of alternative pricing methods}\label{sec:littReview}
As explained in the introductory section above, the numerical methods for option pricing fall into the following two categories: differential methods and integral methods. Differential methods provide the solution to the pricing problem by numerically approximating the associated partial differential equation (see the references of the previous section). Integral methods approximate the conditional expectation, giving the arbitrage-free value of a financial contract. Some methods lie at the boundary of the two groups.
The binomial tree technology, and lattice methods in general (\cite{cox1979option}, \cite{Broadie:1996aa}, \cite{vellekoop2006efficient}),
although similar to finite differences in their implementation (\cite{rubinstein2000relation} shows that they are basically the same), belong to the integral family of methods, given that they provide a discrete approximation of a conditional expectation.
\cite{Barone-Adesi:1987aa} start from the partial differential equation representation in the Black-Scholes case, and contribute a closed-form approximation for an American option with continuous dividend yield.

In this section, we focus on integral methods because the main recent contributions, including ours, fall within this family. The integral representation arises from the observation that option valuation in arbitrage-free economies amounts to using linear operators that assign prices today to payoffs at futures dates. These linear operators are conditional expectation operators. In multi-period economies, the temporal consistency in valuation ensures that the family of such operators satisfies a semigroup property (\cite{Garman:1985aa}). Because most valuation models correspond to Markov environments for the price of the underlying assets, the semigroup property is formalized by the law of iterated expectations restricted to Markov processes (\cite{Hansen:2009aa}), where the Markov states are indexed by the time horizons of potential exercise. From a computational perspective, pricing algorithms attempt to approximate the following: (i) the conditional expectation with respect to the information at one potential exercise date and (ii) the recursive relationship between conditional expectations at different potential exercise dates.
\cite{Geske:1984aa} express the problem of pricing an option with discrete exercise dates as a multivariate integral.
\cite{Kim:1990aa}, \cite{jamshidian1992analysis}, and \cite{carr1992alternative} extend this approach to continuous time and express the early exercise premium of American contracts as an integral term whose interpretation is the risk-neutral valuation of a continuous cash flow. This additional integral is a function of the exercise boundary.
\cite{Bunch:2000aa}, \cite{huang1996pricing}, \cite{Ju:1998aa}, and \cite{ibanez2003robust} contribute algorithms to compute the free boundary and the option prices in the Black-Scholes model environment.
\cite{detemple2002valuation} extend these contributions to diffusion processes with stochastic interest rates. These approaches cannot easily accommodate discrete dividends.
\cite{roll1977analytic}, \cite{geske1979note}, and \cite{whaley1981valuation} provide a closed-form approximation of pricing an
American option paying a single discrete dividend in the Black-Scholes framework. \cite{medvedev2010pricing} develop an approximation of American option prices under stochastic volatility and stochastic interest rates using short-term asymptotics (for the Black-Scholes case, see \cite{Lamberton:2003aa}). These approximations are accurate and flexible, but cannot accommodate discrete dividends. They are mainly suited for index and exchange rate options.

Another approach in approximating the pricing operator is simulating the trajectories of the underlying asset to compute the conditional expectation. One of the first simulation-based analyses of American options in the Black-Scholes case is that of \cite{broadie1997pricing}. A recent stream of Monte Carlo methods started with the work of \cite{longstaff2001valuing}, developing into the so-called duality approach of \cite{haugh2004pricing}, \cite{rogers2002monte}, and \cite{andersen2004primal}. The duality approach method uses simulations to find a lower bound and an upper bound of the true price.
In this sense, it is an extension of \cite{Broadie:1996aa}. As shown in the numerical experiments of \cite{andersen2004primal}, the interval can be tight, which makes the method very precise. The advantage of this approach is that it can handle virtually any type of process dynamics, state variable structure, and payoff specification. However, to achieve precision, as in all simulation-based methods, the discretization step in the time dimension must be small (except trivial cases like the geometric brownian motion), making computations lengthy. \cite{desai2012pathwise} refine the duality approach by contributing a new algorithm to compute the upper bound, improving it in terms of speed.

We can express conditional expectations in terms of series developments on appropriate basis functions, which is the approach we follow when projecting the Green and value functions on localized basis functions in the subsequent sections. We can already find simple and tractable projections of complex payoffs in the literature via the representation on a set of basis functions (such as polynomials in \cite{madan1994contingent}, \cite{lacoste1996wiener}, \cite{darolles2000approximating}).
\cite{chiarella1999evaluation} suggest a fast recursion projection method based on Hermite polynomials in the Black-Scholes model (for an extension of the model in \cite{merton1976option}, see \cite{chiarella2005pricing}).
They also explain how their method provides a viable numerical method for the implementation of the path integral approach to option pricing as described in \cite{linetsky1997path}, for example. Our paper is an extension of the former in that we use sampling instead of projections to improve speed. This family of models typically assumes a geometric Brownian motion for the underlying asset, with some exceptions accommodating for jumps but not for stochastic volatility.

In a more recent application of Hermite polynomials to pricing,  \cite{xiu2014hermite} provides a general, closed-form approximation of European, but not path-dependent, option prices by using a finite-term expansion of the transition density (for related work on expansions, see also \cite{ait2007maximum, ait2010estimating}). \cite{kristensen2011adding} approximate the option price by expanding the difference between the true model price and the Black-Scholes price.
They can price some specific path-dependent options, such as barrier options but cannot address more general path-dependent payoff features, such as discrete dividends.
Our recursion projection method shares some features with the dynamic programming approach of \cite{ben2009dynamic} in that they also approximate the value function on a fixed grid and interpolate with local polynomials to reconstruct the value surface. They assume a GARCH process for the stochastic volatility and obtain closed-form formulas for the conditional expectations to compute the grid of values at a previous time. Though computationally efficient, their method is restricted to the GARCH family of processes and does not naturally accommodate dividends.

As explained above, our method shares features with the quadrature family of models. In addition to the contributions cited above, \cite{jackson2008fourier} introduce a Fourier space time stepping algorithm to price Bermudan options under Levy processes. This efficient method is not extendible to the stochastic volatility framework.

\section{An introductory example}\label{sec:example}

In this section, we show how a pricing problem for a European payoff in the Black-Scholes model translates into a functional projection. The derivation uses elementary calculus concepts. Then, we explain how we exploit the projection to build fast recursive schemes to value path-dependent options such as Bermudan and American options when the stock pays discrete dividends. We design this introductory example to emphasize the intuition underpinning our approach.

\subsection{Description of the method: the European case }\label{sec:example_european}
Let $S_t$ be the price of the underlying asset at date $t$ and assume that interest rates are constant to facilitate
exposition. For $S_t = x$, the value function $V(x,t)$ for a European option is given by the following
conditional expectation:
\begin{equation}\label{eq:condexp}
	V(x, t)=\mathbb{E}\big[e^{-r(T-t)}H(S_T, T)\big|S_t = x\big],
\end{equation}
where $H(S_T,T)$ is the payoff function expressed as a function of time $T$ and of the value of
the underlying asset $S_T$ at maturity date  $T$, and $r$ is the constant risk-free interest rate. When the  pricing operator in \eqref{eq:condexp}
admits a state price density $G(x,t; y, T)$, the so-called Green function, which is the discounted value of the transition probability density from point $x$ at time $t$
to point $y$ at time $T$, we obtain the familiar integral form:
\begin{equation}\label{eq:integral}
	V(x, t) = \int \mspace{-3mu}  G(x,t; y, T) H(y, T) dy.
\end{equation}
Now consider a regularly spaced grid of points $\{y_1, y_2, \dots, y_N\}$ that defines a finite interval $D=[y_1-\Delta y/2, y_N +\Delta y/2]$, with $\Delta y= y_{j} - y_{j-1}$.
We know that, under appropriate regularity conditions, the integral \eqref{eq:integral} restricted to the interval $D$ can be approximated by the Riemann sum as follows:
\begin{equation}\label{eq:Riem_sum}
	V(x, t) \sim \sum_{j=1}^{N-1}  G(x, t; y_j, T)\, H(y_j, T)\,\Delta y, \qquad t<T,
\end{equation}
where the $\sim$ symbol means that the right hand term converges to the left hand term as \mbox{$\Delta y \to 0$}.
The representation \eqref{eq:Riem_sum} is known as the `rectangle method' in standard integral calculus.
If we are interested in computing the value of $V(x,t)$ on a regularly spaced grid of points
$\{x_1, \dots, x_M\}$, for instance to plot the value function of the contract,
or to compute the Greeks, we can express \eqref{eq:Riem_sum} in a matrix form as follows:
\begin{equation}\label{eq:matrix}
	\mathbf{V}(t)  \sim  \mathbf{v}(t)=\mathbf{G}(t;T)\mathbf{H}(T),
\end{equation}
where \mbox{$\mathbf{V}(t) =  \big(V(x_1,t), \dots, V(x_M,t))'$},  $\mathbf{v}(t)$ is the approximation of $\mathbf{V}(t)$  obtained from the $N \times 1$ vector
$\mathbf{H}(T)$  with entries
 $\mathbf H_j= H(y_j,T)$
and the $M\times N$ matrix $\mathbf{G}$  with entries \mbox{$\mathbf{G}_{ij} = \Delta y\, G(x_i,t; y_j,T)$}.
Equation \eqref{eq:matrix} describes a discretization of the economy.
We can think of the matrix $\mathbf{G}$ as a matrix of Arrow-Debreu prices,
in which the rows represent discrete states $\{x_i\}_{i=1,\dots,M}$ at
the current date $t$, and the columns represent discrete states
$\{y_j\}_{j=1,\dots, N}$ at the future date $T$.
We can interpret the column vectors $\mathbf{V}(t)$ and $\mathbf{H}(T)$ as vectors of state contingent values at dates $t$ and $T$, respectively.
The matrix operator $\mathbf{G}$ discounts future payoffs at date $T$ to current prices at date $t$ (\cite{Garman:1985aa}).

Let us more closely examine the elements of $\mathbf{H}(T)$. For readability, we define
$\underline{y}_i = y_i - \frac{\Delta y}{2}$ and $\overline{y}_i = y_i + \frac{\Delta y}{2}$.
We can interpret the value $ H(y_j, T)$  as an approximation of the quantity
\begin{equation}\label{eq:InnerProduct}
	(1/\Delta y)\int_{\underline{y}_j}^{\overline{y}_{j} }\mspace{-9mu}H(y,T)dy= (1/\sqrt{\Delta y})\int \frac{1}{\sqrt{\Delta y}} \mathbb{I}_{\underline{y}_j, \overline{y}_{j} } H(y,T)dy \stackrel{\rm def}{=} \frac{1}{\sqrt{\Delta y}}\int\mspace{-3mu} e_{j}(y)H(y,T)dy,
\end{equation}
where $e_{j}(y)=\frac{1}{\sqrt{\Delta y}}\mathbb{I}_{\underline{y}_j, \overline{y}_{j} }$, and where $\mathbb{I}_{\underline{y}_j, \overline{y}_{j} }$ is the indicator function of the interval $[\underline{y}_j, \overline{y}_{j})$.  Because
$\{e_{j}(y)\}_{j=1,\dots,N}$ is an orthonormal set given
the standard $L^2$ inner product $\langle f,g\rangle = \int \mspace{-6mu}f(x)g(x)dx$, we can view the entries of the vector $\mathbf{H}(T)$ as an approximation
of $(1/\sqrt{\Delta y})$ times the coefficients of the decomposition of the payoff function $H(y, T)$
on the set of orthonormal indicator functions defined by the grid $\{y_1, \dots,y_N\}$.
A similar argument can be applied to the coefficients of the $\mathbf{G}(t;  T)$ matrix: every
row of $\mathbf{G}(t;T)$ is given by an approximation of $\sqrt{\Delta y}$ times the coefficients of the projection of the conditional density
$G(x_i, t; y, T)$ on the same orthonormal set.
The different choices in  the normalization factors for the entries of $\mathbf{G}$ and $\mathbf{H}$
are justified by our wish to interpret all of the quantities appearing in \eqref{eq:matrix} as prices.
We have that $V(x,t)$ and $H(y,T)$ are value functions and $G(x,t;y,T)$ is a state price density.

Overall, we can interpret the numerical approximation of the integral
 \eqref{eq:integral} as a projection of the functions on an orthonormal basis.
Such an interpretation is key
in the generalization of the recursive projection approach to more sophisticated models in the next section. In our case, the functional projection boils down to sampling the given functions
on a grid of $N$ points $\{y_j\}_{j=1,\dots,N}$ via a discrete transform.
From a computational perspective, the entries of $\mathbf{H}(T)$ summarize how payoffs depend on the price of the underlying assets at future exercise dates,
and the entries of $\mathbf{G}(t;T)$ summarize how the price of the underlying assets transits from one state to another according to the elapsed time between exercise dates. Figure \ref{fig:projection-1d} gives a graphical representation of the two computational steps of our fast recursive projection approach, (i) the projection step, and (ii) the recursive step. The value function at $t$, on the right,
 and the state price density for a given value of $x=S_{t}$, on the left, are sampled (the projection step); the obtained arrays of values are multiplied element by element, and the products are summed to obtain the value of
 $V(x, t)$ (the recursive step).

	\medskip
	{\centering
	[Figure  \ref{fig:projection-1d} about here]\medskip

	}

\subsection{Description of the method: the path-dependent case }\label{sec:exampl_PathDep}
Let us now address the valuation of path-dependent contracts. We start by considering a Bermudan option. We consider a set $\{t_1=t, \dots, t_L=T\}$ of exercise dates. At each $t_l$, the holder of a Bermudan option may decide whether to exercise. He exercises if the intrinsic value \mbox{$H(S_{t_l}, t_l) \stackrel{\rm def}{=} (S_{t_l} - K)_+= \max\{S_{t_l} -K,0\}$} is higher than the value of keeping the option, i.e., the continuation value. Bermudan options are the ideal building blocks for studying American call
options on dividend-paying stocks. It is well known that it can be optimal to exercise American call options immediately before ex-dividend dates
$\{t_h, t < t_h < T\}_{h=1,\dots, H}$; for instance, see \cite{pool2008failure} for a discussion on early exercise strategies. The implication is that we must monitor the option value function $V(x,t)$ immediately before the ex-dividend dates, when the intrinsic value $(S_{t_h-\epsilon} - K)_+$  for a small $\epsilon >0$ can be larger than the continuation value $V(S_{t_h},t_h)$. Then an American call option shares with a Bermudan option the feature that its value function must be evaluated only at a finite number of dates.

 The semigroup property of the pricing operator ensures that we compute the value function $V(x,t)$ of a Bermudan option recursively.
The recursion consists of moving backwards in time and computing at each $t_l, l=1, \dots , L-1$:
\begin{equation}\label{eq:rec_exp}
	V(x,t_l) = \max\bigl\{H(x,t_l), \mathbb{E}\bigl[e^{-r(t_{l+1} - t_{l})}V(S_{t_{l+1}}, t_{l+1})|S_{t_l}=x\bigr]\bigr\},
\end{equation}
with the boundary condition $V(y, t_L)=H(y,t_L) $.
To speed up the recursion, we impose the condition that the grid of values  $\{y_j\}_{j=1,\dots,N}$ at which we sample the function $V(y, t_{l+1})$  and the grid  $\{x_i\}_{i=1,\dots, M}$ at which we compute the function $V(x, t_l)$, coincide at each exercise date, which means that  $M=N$. From now on, we use the $x$ variable as a generic conditioning value, e. g. the value of the underlying at date $t$, as in $V(x,t)$. If $x$ takes a specific value in the grid $\{y_{i}\}_{i=1,\dots,N}$, then we write $V(y_{i},t)$.
Then, in the matrix notation of the approximation, we obtain the following:
\begin{equation}\label{eq:rec_matr}
	\mathbf{V}(t_l)  \sim \mathbf{v}(t_l)= \max \bigl\{\mathbf{H}(t_l), \mathbf{G}(t_l;t_{l+1})\mathbf{v}(t_{l+1})\bigr\},
\end{equation}
and the approximation $\mathbf{v}({t_l})$ is the input for the following recursion step to compute an approximation $\mathbf{v}({t_{l-1}})$  of  $\mathbf{V}({t_{l-1}})$.
Figure \ref{fig:projection-1dmatrix} shows how combining projection and recursion steps with identical grids at each exercise date translates the pricing problem of path-dependent options into a sequence of matrix time vector multiplications.

	\medskip
	{\centering
	[Figure  \ref{fig:projection-1dmatrix} about here]\medskip

	}

The convergence properties of $\mathbf{v}({t_l})$ to $\mathbf{V}({t_l})$  are formally established in Proposition \ref{prop:convBS} of Section \ref{sec:RecProj}.
From \eqref{eq:rec_matr}, it is clear that computations only occur at the exercise dates defined in the Bermudan contract, and do not require any input at any other point in time.
Furthermore, under the previous additional constraints on the grids, if the time
interval $\tau = t_{l+1} - t_l, l=1, \dots, L-1$, is constant and if the pricing operator
enjoys a stationarity property (time translation invariance), then the matrix $\mathbf{G}(t_l;t_{l+1})=\mathbf{G}(\tau)$
has constant entries, and the algorithm only involves one single computation of the matrix.

The methodology easily extends  in the presence of discrete dividends
paid on potential exercise dates\footnote{We can easily extend to the case in which the ex-dividend date does not belong to the set of exercise dates. We do not explicitly state this case, given that we are mainly interested in Bermudan options, as a building block for studying American options.}.
The implication is that the set of ex-dividend dates is a subset of the Bermudan exercise dates and that we have
$\{t_{h}\}_{h=1,\dots,H} \subset \{t_{l}\}_{l=1,\dots, L}$.
    We only need to add the dividend $\delta$ to the continuation value in  (\ref{eq:rec_exp}).
Hence, in order to price an American option on a dividend-paying stock, Equation \eqref{eq:rec_matr} must be modified by sampling the state price density $G(x,t_{l};y,t_{l+1})$ at
the grid $\{y_i-\delta(y_{i})\}_{ i=1,\dots, N}$ for the conditioning  value $x$, whenever $t_{l}\in \{t_{h}\}_{h=1,\dots,H}$. The entries of the matrix $\mathbf{G}(t_{l}; t_{l+1})$ then become $\mathbf{G}_{ij}=G(y_i-\delta, t_l; y_j, t_{l+1})\Delta y$. Given the freedom in choosing where to sample $G$, $\delta(x)$ could be any function of $x$. If $\delta(x)=r_d \, x$, then we can accommodate for a proportional dividend. If $\delta(x)=d$, then we
can accommodate for  a discrete dividend amount $d$\ \footnote{When we consider discrete dividends, we rule out arbitrage situations where the dividend is too large with respect to the stock price (see \cite{haug2003back} for a discussion).}. If $\delta(x)=0$, then we return to the Bermudan option case. The value function $\mathbf{V}(t_{h})$ still gives the value of the contract at the grid points $\{y_{1}, \dots, y_{N}\}$; thus we can use its approximation $\mathbf{v}(t_{h})$ as the input for the following step of the algorithm, and the recursive property of the algorithm is maintained.
Figure \ref{fig:dividends}, shows how the recursive scheme changes to accommodate for dividends.

	\medskip
	{\centering
	[Figure  \ref{fig:dividends} about here]\medskip

	}

The early exercise decision for American long put holders is more complicated. We consider a Bermudan put with exercise dates $\{t_1,\dots,t_L\}$, and assume that the dates $\{t_{h}\}$ at which dividends are paid form a subset of the exercise dates. Then, by taking $L$ large,
our approach also provides a quick approximation for pricing American put options on single stocks
paying discrete dividends.
In the Black-Scholes model, the state price density is known in the following closed form:
\begin{equation}\label{eq:GreenBS}
	G(x,t; y, T)= \frac{1}{y}\frac{e^{-r(T-t)}}{\sqrt{2\pi \sigma^2(T-t)}}\exp\biggl(-\frac{\bigl(\log y - \log x - (r-\sigma^2/2)(T-t)\bigr)^2}{2 \sigma^2(T-t)}\biggr),
\end{equation}
where $\sigma$  is the volatility and the pricing operator satisfies the stationarity property. Hence, we expect to obtain a fast, simple and accurate numerical algorithm.
Before reporting the numerical results, we stress that we have kept the setup of this introductory example as simple as possible to emphasize intuition by limiting technical details. Although the computational results already speak in favor of our construction, it is in the following sections that the recursive projections method deploys its full potentiality with more complex price dynamics.

\subsection{Numerical illustrations in the Black-Scholes model}\label{sec:BSillustration}

As a first numerical example in the Black-Scholes framework, we compare the convergence speed of a binomial tree and of the recursive projections method in pricing an American call option on a dividend-paying stock\footnote{All of the codes are written in C++. The codes are available from the authors upon request.}.
Two popular modeling choices for the dividend payment are a known  cash amount $d$ or a known dividend yield $r_d$. The latter is computationally friendly because it leads to a recombining tree. However, the empirical evidence shows that corporations tend to commit to paying out fixed amounts at regular dates and to smooth their dividends rather than adjusting them  downwards and signaling a decrease in cash flows (for a signaling based theory on dividend policy see, for instance, \cite{Miller1985zr}). These arguments justify our preference for modeling dividends as fixed known amounts rather than as given yields. The known dividend amount assumption does not lead to a recombining tree, and a new tree is originated at each node following an ex-dividend date, increasing the numerical complexity of the problem. The work of \cite{vellekoop2006efficient} provides a recent enhancement of the classical binomial tree method which incorporates discrete dividend payments through an approximation of the continuation value of the option at the ex-dividend dates. This new algorithm has been proven to be substantially faster than the standard non-recombining binomial tree, and is therefore a reliable benchmark for this simulation exercise.

	\medskip
	{\centering
	[Figures  \ref{fig:Vellekoop} and  \ref{Figure:BS_amcall_div} about here]\medskip

	}

Figure  \ref{fig:Vellekoop} compares the convergence speed of the enhanced binomial tree and that of the recursive projections method in pricing an American call option on a discrete dividend-paying stock. The option has a maturity of $T=3$ years and a dividend $d=2$ is paid out at the end of each year.  Other parameters, namely the interest rate, volatility and strike price,  are set equal to $r=0.05$, $\sigma=0.2$,  and $K=100$, respectively. We compute 3 prices: at-the-money, in-the-money and out-of-the-money, corresponding to $S_0=80, 100$, and $120$, respectively.
The \emph{true} values of 7.180, 18.526, and 34.033 are obtained with $10000$ time steps in the binomial tree. The graphs show that, across the three different values of $S_{0}$, the recursive projections enjoy an increase of speed of approximately a factor $10$ for a comparable level of precision. The speed advantage is even larger if we consider that a new tree is needed for each value of $S_{0}$. Instead, the recursive projections method delivers the entire value function $\mathbf{v}(0)$ at once in a straightforward manner. This feature is particularly useful in computing Greeks through numerical differentiation.
As an additional benchmark, Figure \ref{Figure:BS_amcall_div} displays the convergence speed of the recursive projections jointly with the one of a standard non-recombining tree. Even though the non-recombining tree is known to be an inefficient method, it is still used as a common reference point in the literature, and we show this graph for comparison purposes. We can see that the gain of speed of the recursive projection is of the order of $10^{4}$.
In Section E of the supplementary online Appendix, we extend the above pricing exercise, and further compare the recursive projections with the Monte Carlo based method by \cite{longstaff2001valuing}. The recursive projections are faster than the Longstaff Schwartz method by at least four orders of magnitude.
As an aside,  for $S_{0}=100$, if we approximate the known constant dividend $d=2$ with a known continuous dividend yield\footnote{The yield is obtained by considering the dividends paid at $t=1$ and $t=2$ only, because the dividend paid at $t=3$ has no impact on the price of the option. Considering a dividend yield of $2\%$ would provide an option value of 16.857, which is a much larger error.} $r_d=0.013$, then a binomial tree with $10000$ steps delivers a value of 18.213 instead of 18.526, with a relative error of approximately $169bp$.
This error is far above observed bid-ask spreads. This simple example points to the importance of using models that can explicitly address discrete dividends in empirical analysis, instead of using approximations based on continuous dividend yields.
In Section \ref{sec:boundaries_empirics}, we more extensively analyze the impact on the exercise boundary of the choice of a continuous dividend yield or a discrete dividend.

Moreover, we have chosen a sampling scheme that is equivalent to projecting the payoff function on a set of basis functions that are well localized, in the sense that their support is a closed interval. The implication is that local features of the payoff function, such as a discontinuity, are described by the coefficients relative to one or at most two basis functions, those lying next to the discontinuity.
This description avoids a noisy approximation induced by spurious oscillations when projecting discontinuities on basis functions defined on the entire domain, such as the Fourier sine-cosine basis or the Hermite polynomial basis.
From a computational perspective, this property translates into an accurate approximation even for payoffs with strong discontinuities, such as a digital payoff  $H(S_{t_l}, t_l) = \mathbb{I}_{S_{t_l>K}}$ in a Bermudan digital call option. The discontinuity may introduce noise at most in the coefficient relative to the indicator function of the interval in which the discontinuity is located.
The noise is completely eliminated if we make sure that that the strike value lies in between two consecutive grid points, so that the discontinuity is at the boundary between two consecutive indicator functions.
In this numerical example, we use the standard binomial tree as a benchmark, since the method of \cite{vellekoop2006efficient} provides no advantage in the absence of dividends.
Figure \ref{Figure:BS_berm_digitalcall} (see the caption of the table for the values of the parameters of the example) shows that the binomial tree has problems capturing the discontinuity in the payoff function. Consequently, an extremely slow convergence of the tree method for at-the-money Bermudan digital call options is yielded. The recursive projections are also at least an order of magnitude faster in pricing the out-of-the-money options. The apparent non-monotonic convergence of the binomial tree for $S_{0} = 120$ is because both methods achieve a quick convergence for in-the-money options, and the graph only displays small oscillations on the order of half a basis point around the true value.

	\medskip
	{\centering
	[Figure \ref{Figure:BS_berm_digitalcall} about here]\medskip

	}

\section{Valuation by fast recursive projections}\label{sec:RecProj}

In this section, we generalize the approach developed in the introductory example of Section \ref{sec:example_european}. We consider a set of exercise dates $\{t_l\}_{l=1,\dots, L}$, with $t_{L}=T$.
The starting point is the following identity:
\begin{align}\notag
	V(x,t)
&=
 \int \mspace{-3mu} G(x,t ; y, T)H(y,T) dy
=\mspace{-9mu}
\int \mspace{-6mu}
\Bigl(\sum_{i=-\infty}^{\infty} \mspace{-6mu}
\langle G , \varphi_{i}\rangle  \varphi_{i}(y)\Bigr)
\Bigl(\sum_{j=-\infty}^{\infty} \mspace{-6mu}
\langle H , \varphi_{j}\rangle \varphi_{j}(y)\Bigr)
 dy \\
&=\mspace{-9mu}
\sum_{i,j=-\infty}^{\infty} \mspace{-6mu}
\langle G , \varphi_{i}\rangle
\langle H , \varphi_{j}\rangle
 \mspace{-6mu}
 \int \mspace{-6mu} \varphi_{i}(y)\varphi_{j}(y) dy  \label{eq:generalprojection}
	 =
 \sum_{j=-\infty}^{\infty}
\langle G,\varphi_{j}\rangle
\langle H ,\varphi_{j}\rangle,
\end{align}
where $\{\varphi_{j}\}_{i\in \Z}$ is a generic localized orthonormal basis in $L^2$,  that is the support of each function $\varphi_{j}$ is a finite closed interval.  $\langle G(x,t ; \cdot, T) \varphi_{j}(\cdot)\rangle$ and $\langle H(\cdot, T),\varphi_{j}(\cdot)\rangle$ (abbreviated in $\langle G,\varphi_{j}\rangle$ and $\langle H, \phi_{j}\rangle$ respectively) are the coefficients of the linear projection of $G$ and $H$ on $\{\varphi_{j}\}_{i\in \Z}$. The valuation by fast recursive projections is based on two steps: a projection step and a recursive step.
The projection step consists of projecting $G$ and $H$ at time $T$ on the basis $\{\varphi_{j}\}_{j\in \Z}$ and computing the coefficients $\langle G,\varphi_{j}\rangle$ and $\langle H,\varphi_{j}\rangle$. Due to finite support of the localized basis functions $\{\varphi_{j}\}_{j\in \Z}$, all the coefficients $\langle H,\varphi_{j}\rangle$ computed as the inner product between the function $H$ and the localized basis function $\varphi_{j}$ are automatically well defined, even if $H$ is not $L^2$-measurable. The recursive step consists of transmitting the coefficients back in time by computing the final sum in Equation (\ref{eq:generalprojection}).
We develop the method by considering the orthonormal set of indicator functions $\{e_{j}\}_{j\in \Z}$, $e_{j}=1/\sqrt{\Delta y}\, \mathbb{I}_{\underline{y}_{j}, \overline{y}_{j}}$ for an equally spaced grid $\{y_{i}\}$ with step $\Delta y$ and $\underline{y}_{j}=y_{j}-\Delta y/2$, $\overline{y}_{j}=y_{j}+\Delta y/2$ . The set $\{e_{j}\}_{j\in \Z}$ is a basis in $L^2$ when $\Delta y \to 0$.
Then, as shown in Equation \eqref{eq:InnerProduct}, we can approximate the coefficients $\langle H, e_{j}\rangle$ and $\langle G, e_{j}\rangle$ with $\sqrt{\Delta y}H(y_{j},T) $ and  with $\sqrt{\Delta y}\, G(x,t ; y_{j}, T)$. Thus, the projection step simply consists of sampling the relevant functions at the grid points $\{y_{j}\}_{j\in \Z}$.
The recursion step stems from the observation that, if $x$ in Equation \eqref{eq:generalprojection} takes values in the grid $\{y_{j}\}_{j\in \Z}$, then $v(y_{i},t_{L-1}) = \sum_{j=-\infty}^{\infty} G(y_{i},t_{L-1} ; y_{j}, T)H(y_{j}, T)\Delta y$ and $v(y_{i},t_{L-2}) = \sum_{j=-\infty}^{\infty} G(y_{i},t_{L-2} ; y_{j}, t_{L-1})v(y_{j}, t_{L-1})\Delta y$ are approximations of $V(y_{i}, t_{L-1})$ and $V(y_{i}, t_{L-2})$. That is, there is a recursive linear expression that connects the value of the contract at consecutive times and at different points of $\{y_{j}\}_{j\in \Z}$.
Here, all the functions are sampled at the same grid $\{y_{j}\}_{j\in \Z}$ for each potential exercise date $\{t_{l}\}_{l=1,\dots,L}$. Doing so allows the matrix representations of Equation \eqref{eq:rec_matr}.
In the following sections, we characterize the projection step and the recursive step in two different frameworks.
In Sections \ref{sec:ProjectionBS} and \ref{sec:RecursBS}, we consider the case in which the state price density is known in closed form, as in the Black-Scholes example developed in Section \ref{sec:example} or in the Merton jump-diffusion model. In Sections \ref{sec:ProjectionHes} and \ref{sec:RecursHes}, we characterize the projection step and the recursion step when the characteristic function of the state price density has a known analytic form, as in the Heston model.
In our presentation, we generically refer to this second class of models as to the ``stochastic volatility'' case. The methodology developed in Sections \ref{sec:ProjectionHes} and \ref{sec:RecursHes} covers affine jump-diffusion models and Levy models, in which we price by transform analysis.
Here, the key requirement is to be able to numerically compute the characteristic function of the state price density. In Section \ref{sec:NumericalHeston}, we provide numerical examples in the Heston case.

\subsection{Merton-Black-Scholes: the projection step}\label{sec:ProjectionBS}

The projection step is based on an approximation of the payoff function and of the state price density by the set of orthonormal functions $\{ e_{j} (y) \}_{j\in \Z}$:
\begin{align} \label{eq:projection}
	\tilde H(y,T) &= \sqrt{\Delta y}\, \sum_{j=-\infty}^{\infty} H(y_{j},T) e_{j}(y), \\ \label{eq:projection2}
	\tilde G(y_{i},t; y, T) & = \sqrt{\Delta y}\, \sum_{j=-\infty}^{\infty} G(y_{i},t; y_{j}, T) e_{j}(y).
\end{align}
The values of $\tilde H(y,T)$ and $\tilde G(y_{i}, t; y, T)$ in each interval $[\underline{y}_{j}, \overline{y}_{j})$ are given by the values of the functions $G$ and $H$ sampled in the mid point $y_{j}$ of each interval. Then, we have that $\tilde H(y,T)$ and $\tilde G(y_{i}, t; y, T)$ are piecewise constant approximations of $H(y,T)$ and $G(y_{i}, t; y, T)$.
In the algorithm, we only need the quantities $\{H(y_{j},T)\}_{j\in \Z}$ and  $\{G(y_{i},t; y_{j}, T)\}_{j\in \Z}$.
We use the full representations \eqref{eq:projection} and \eqref{eq:projection2} in Section A of the supplementary online Appendix to prove Proposition \ref{prop:convBS}.
There, we show that relying on the sampling approximations \eqref{eq:projection} and \eqref{eq:projection2} instead of the orthogonal projections given by the inner products as in Equation \eqref{eq:generalprojection}, leaves the convergence properties of the pricing algorithm unaffected.

\subsection{Merton-Black-Scholes: the recursive step}\label{sec:RecursBS}

The following proposition gives a recursion formula that relates the approximated values of the option at different points of the grid $\{y_{j}\}_{j=-\infty}^{\infty}$ and at different points in the set of exercise dates $\{t_{l}\}_{l=1,\dots,L}$. It also states the rate of convergence of the algorithm.

\begin{proposition}\label{prop:convBS}
Let $H(y,T)$ be such that $|H(y,T)-H(y',T)|<C|y-y'|$ for a positive constant $C$, and for $|y-y'|<\Delta y$.
Furthermore, let $v_{i}(t_{l})$ be defined for a set of dates $\{t_{l}\}_{l=1,\dots,L}$, with $t_{L}=T$, as follows:
\begin{alignat}{3}\label{eq:recursiveBS}
	v_i (t_{l}) &= \max\Bigl\{H(y_{i}, t_{l}), \sum_{j \in \Z} G(y_{i}, t_{l};y_{j}, t_{l+1})H(y_{j}, t_{l+1}) \Delta y \Bigr\},		&	\quad  	& \mathrm{for}	& \quad l&=L-1 ,	\\
	\label{eq:recursiveBS2}
	v_i (t_{l}) &= \max\Bigl\{H(y_{i},t_{l}), \sum_{j\in \Z} G(y_{i}, t_{l};y_{j}, t_{l+1})v_{j}(t_{l+1})\Delta y \Bigr\}, 	&			& \mathrm{for}	& l&=1,\dots, L-2.
\end{alignat}
Then, for each $t_{l}$ in $\{t_{1},\dots, t_{L-1}\},$ the approximated values $v_{i}(t_{l})$ defined in \eqref{eq:recursiveBS} and \eqref{eq:recursiveBS2} converge to the true value $V(y_{i}, t_{l})$ with an approximation error of the order $O\bigl((\Delta y)^{2}\bigr)$.
\end{proposition}
\begin{proof}
See Section A of the supplementary online Appendix.
\end{proof}

In principle, the continuity condition on $H(y,T)$ rules out digital payoffs. We need the condition to hold within each interval $\Delta y$ centred on the grid $\{y_{j}\}_{j=-\infty}^{\infty}$. We can still price digital options by ensuring that the strike value lies in between two consecutive grid points. Then the convergence properties stated in Proposition \ref{prop:convBS} remain true.
This procedure is not the same as placing nodes in the quadrature method because the grids remain the same for all dates  $\{t_{l}\}_{l=1,\dots,L}$
provided that the strike price does not change with time.
With a fixed grid, we can price  even more exotic options, e.g. a digital call with a down-and-out feature, provided that the barrier occurs at the same value of the underlying for each $t_{l}$.

The main difference between Equations \eqref{eq:recursiveBS} and \eqref{eq:recursiveBS2} is provided by the following explanation. In the right-hand side of  \eqref{eq:recursiveBS}, we find the exact values taken by the payoff function $H(y,T)$ on the grid $\{y_{j}\}_{j\in \mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}}\newcommand{\E}{\mathbb{E}}$, and there is no approximation of the payoff. On the right-hand side of Equation \eqref{eq:recursiveBS2}, we find the values  $\{v_{j}(t_{l+1})\}_{j\in \mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}}\newcommand{\E}{\mathbb{E}}$ obtained in the previous step of the algorithm, and these are approximations of the true values $\{V(y_{j}, t_{l+1})\}_{j\in \mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}}\newcommand{\E}{\mathbb{E}}$.
Regardless of this fundamental difference, Proposition \ref{prop:convBS} states that the convergence rate is the same for both cases.
In the case of a European option with maturity $T$, by taking $t=t_{l}$ in \eqref{eq:recursiveBS}, we obtain the approximated price of the option at $t$. Equation  \eqref{eq:recursiveBS2} allows us to recursively compute the values of the option at different points in time, and thus to price Bermudan options, American options, and other types of path-dependent options.

In the implementation of Section \ref{sec:example_european}, we project the payoff function on a finite number of indicator functions. Thus, the values $\{H(y_{j}, T)\}_{j=1,\dots,N}$ and $\{\Delta y\, G(y_{i},t; y_{j}, t_{l+1})\}_{j=1,\dots,N}$ are then the entries of the finite dimensions matrices $\mathbf{H}(T)$ and $\mathbf{G}(t_{l};t_{l+1})$, respectively. Hence, we can express Equations \eqref{eq:recursiveBS} and \eqref{eq:recursiveBS2} in exactly the same matrix forms as in \eqref{eq:matrix} and \eqref{eq:rec_matr}.
By choosing $N$ sufficiently large, we can make the error introduced by truncating the infinite summations in  \eqref{eq:projection} and \eqref{eq:projection2} arbitrarily small. The truncation only suppresses the indicator functions localized in regions where the state price density vanishes, and where we can neglect the local contribution to the computed expectation\footnote{Truncating $N$ is equivalent to choosing an $y_{min}$ and $y_{max}$, and to restrict our analysis to $L^2([y_{min}, y_{max}])$.
This is the standard assumption in each quadrature based approach. This also insures that $H\in L^2([y_{min}, y_{max}])$, even though the proofs in Appendices A and B of the supplemental file do not need the measurability of $H$.}.
When the time step is a constant $\tau=t_{l+1} - t_{l}$, we obtain the fast and easily implementable algorithm of the introductory example.

Furthermore, as in Section \ref{sec:exampl_PathDep}, whenever an exercise date $t_{l}$ coincides with a dividend-paying date $t_{h}$,
we only need to replace the entries $\{\Delta y\, G(y_{i},t_{l}; y_{j}, t_{l+1})\}_{j=1,\dots,N}$ of  $\mathbf{G}( t_{l}; t_{l+1})$ with the values $\{\Delta y\, G(y_{i}-d, t_{l}; y_{j}, t_{l+1})\}_{j=1,\dots,N}$, whenever $t_{l}\in \{t_{h}\}_{h=1,\dots,H}$, to accommodate for a discrete dividend $d$.

\subsection{Stochastic volatility: the projection step}\label{sec:ProjectionHes}
In the class of  stochastic volatility models, there are two state variables, the underlying asset $S_{t}$ and the variance $\sigma^{2}_{t}$.
The bivariate state price density $G_{2}(S_{t}, \sigma^{2}_{t}, t; y, w, T)$ describes the discounted transition probability density from the asset level $S_{t}$ and variance level $\sigma^{2}_{t}$ at time $t$ to the asset level $y=S_{T}$ and variance level $w=\sigma^{2}_{T}$ at time $T$.
Its Fourier transform is denoted by $\hat{G}_{2}(S_{t}, \sigma^{2}_{t}, t; \lambda, \kappa, T)$, so that $G_{2}(S_{t}, \sigma^{2}_{t}, t; y, w, T) =\displaystyle \frac{1}{4\pi^{2}}{\iint} d\lambda d\kappa e^{-\iota( \lambda y +\kappa w)}\hat{G}_{2}(x, \xi, t; \lambda, \kappa, T)$, where $\iota$ is the imaginary unit.

For the values taken by the underlying asset, let the grid $\{y_{j}\}_{j\in \mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}}\newcommand{\E}{\mathbb{E}}$ and the orthonormal set $\{e_{j}(y)\}_{j\in \Z}$ be defined as in the Black-Scholes case.  The Fourier transforms of $\{e_{j}(y)\}_{j\in \Z}$ are denoted  by $\{\hat{e}_j(\lambda)\}_{j\in \Z}$ , such that (see Section C of the supplementary online Appendix for the analytic form of $\hat{e}_j(y)$): $e_{j}(y)=\frac{1}{2\pi}\int d\lambda e^{-\iota \lambda y}\hat{e}_{j}(\lambda)$.

For the variance, we use the same equally spaced grid $\{w_{q}\}_{q\in \mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}}\newcommand{\E}{\mathbb{E}}$ for the values taken by both the variables $\sigma^{2}_{t}$ and $w=\sigma^{2}_{T}$, such that $\Delta w= w_{q+1} - w_{q}$. Let $\{\varepsilon_{q}(w)\}_{q\in \mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}}\newcommand{\E}{\mathbb{E}}$ be the normalized indicator functions centered on the grid $\{w_{q}\}_{q\in \mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}}\newcommand{\E}{\mathbb{E}}$ and of width $\Delta w$, and $\{\hat \varepsilon_{q}(\kappa)\}_{q\in \mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}}\newcommand{\E}{\mathbb{E}}$ be their Fourier transforms. Furthermore, let $\{\lambda_{r}\}_{r\in \Z}$ and $\{\kappa_{z}\}_{z\in \Z}$ be two regularly spaced grids of values taken by the transformed variables $\lambda$ and $\kappa$, with constant widths $\Delta \lambda$ and $\Delta \kappa$, respectively.

The approximation of the payoff function $H(y,T)$ is the same as in Equation \eqref{eq:projection}. Because we do not have the analytical form of $G_{2}(S_{t}, \sigma^{2}_{t}, t; y, w, T)$, we cannot directly sample the transition density, thus we must rely on the approximation given below.
The projection step for the bivariate state price density $G_{2}(S_{t}, \sigma^{2}_{t}, t; y, w, T)$ when the conditioning variables take the values $S_{t}=y_{i}$ and $\sigma^{2}_{t} =w_{p}$, is based on the following approximation:
\begin{align}\label{eq:projectionHesG}
	&\tilde G_{2}(y_{i}, w_{p}, t; y, w, T)\\ \notag
	&=\sqrt{\Delta y \Delta w} \mspace{-6mu}\sum_{j=-\infty,q=1}^{\infty} \Bigl( \frac{1}{4\pi^{2}}\mspace{-9mu}\sum_{r,z=-\infty}^{\infty}\mspace{-9mu} \hat{G}_{2}(y_{i},w_{p}, t; \lambda_{r}, \kappa_{z}, T)\hat{e}_{j}(-\lambda_{r})\hat{\varepsilon}_{q}(-\kappa_{z})\Delta \lambda \Delta \kappa\Bigr)e_{j}(y)\varepsilon_{q}(w)	\\ \notag
	& \stackrel{\rm def}{=} \sqrt{\Delta y \Delta w} \mspace{-6mu}\sum_{j=-\infty, q=1}^{\infty}\Gamma_2(y_{i}, w_{p},t; y_{j}, w_{q}, T)e_{j}(y)\varepsilon_{q}(w),
\end{align}
where the second equality in \eqref{eq:projectionHesG} defines the quantities $\{\Gamma_2(y_{i}, w_{p},t; y_{j}, w_{q}, T)\}_{j \in \Z, q \in \mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}}\newcommand{\E}{\mathbb{E}}$.
To parallel the discussion following Equation \eqref{eq:projection2} in the Black-Scholes case, the quantities $\{\Gamma_2(y_{i}, w_{p},t; y_{j}, w_{q}, T)\}_{j \in \Z, q \in \mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}}\newcommand{\E}{\mathbb{E}}$ play the same role as $\{G(y_{i}, t;y_{j}, T)\}_{j \in \Z}$.
We can motivate these approximations as follows. The orthogonal projection of the state price density $G_{2}(y_i,w_{p},t; y,w, T)$ on the two orthonormal sets
$\{e_{j}(y)\}_{j\in \mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}}\newcommand{\E}{\mathbb{E}}$ and $\{\varepsilon_{q}(w)\}_{q\in\mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}}\newcommand{\E}{\mathbb{E}}$ is given by the inner products $\iint \mspace{-3mu}dydw\mspace{1mu} G_{2}(y_i,w_{p},t;y,w,T)e_{j}(y)\varepsilon_{q}(w)$. Because we only know the closed form\footnote{For the closed form of $\hat{G}_{2}(y_i,w_{p}, t; \lambda, \kappa, T)$, see \cite{griebsch2013evaluation}.} of
$\hat{G}_{2}(y_i,w_{p}, t; \lambda, \kappa, T)$ and not $G_{2}(y_i,w_{p},t;y,w,T)$, we exploit the following key relationship:
\begin{align}\notag
	&\iint dydw\mspace{1mu} G_{2}(y_i,w_{p},t;y,w,T)e_{j}(y)\varepsilon_{q}(w) \\ \notag
	& = \iint dydw \frac{1}{4\pi^{2}}{\iint} d\lambda d\kappa e^{-\iota( \lambda y +\kappa w)}\hat{G}_{2}(y_{i}, w_{p}, t; \lambda, \kappa, T)e_{j}(y)\varepsilon_{q}(w) \\ \notag
	&= \frac{1}{4\pi^{2}}{\iint} d\lambda d\kappa\, \hat{G}_{2}(y_i,w_{p},t;\lambda, \kappa, T) \int dy e^{-\iota\lambda y}e_{j}(y)\int dw e^{-\iota\kappa w}\varepsilon_{q}(w) \\
	\label{eq:GammaTwo}
	&= \frac{1}{4\pi^{2}}{\iint} d\lambda d\kappa\, \hat{G}_{2}(y_i,w_{p},t;\lambda, \kappa, T) \hat{e}_{j}(-\lambda) \hat{\varepsilon}_{q}(-\kappa).
\end{align}

\noindent
Each $\Gamma_2(y_{i}, w_{p},t; y_{j}, w_{q}, T)$ is an approximation of the last integral appearing in Equation \eqref{eq:GammaTwo}, obtained by a direct sampling of the Fourier tranforms $\hat G_{2}(x_i,w_{p},t;\lambda, \kappa, T)$, $\hat{e}_{j}(-\lambda)$ and $\hat{\varepsilon}_{q}(-\kappa)$ on the bivariate grid $\{(y_{j},w_{q})\}_{j\in\Z,q\in \mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}}\newcommand{\E}{\mathbb{E}}$, and on the univariate grids $\{\lambda_{r}\}_{r\in \Z}$ and $\{\kappa_{z}\}_{z\in \Z}$.

As in Section \ref{sec:ProjectionBS}, only the quantities $\{\Gamma_2(y_{i}, w_{p},t; y_{j}, w_{q}, T)\}_{j\in \mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}}\newcommand{\E}{\mathbb{E},q \in \Z}$ are the inputs for the pricing algorithm. The representation \eqref{eq:projectionHesG} is only used in Section B of the supplementary online Appendix to prove the convergence properties of the algorithm.

\subsection{Stochastic volatility: the recursive step}\label{sec:RecursHes}
In the stochastic volatility framework, the recursion for a Bermudan option consists of moving backwards in time as in Equation \eqref{eq:rec_exp} with:
\begin{equation}\label{eq:condexpStochVol}
	V(x,\xi, t_{l}) = \max\bigl\{H(x,t_{l}),\E\bigl[e^{-r(t_{l+1}-t_{l})}V(S_{t_{l+1}},\sigma^{2}_{t_{l+1}}, t_{l+1})\big| S_{t_{l}}=x, \sigma^{2}_{t_{l}}=\xi\bigr]\bigr\}.
\end{equation}
Thus, the recursive step in the Heston model is the sampling counterpart of \eqref{eq:condexpStochVol}.
 The following proposition gives a recursion formula that relates the approximated values of the option at different points of the bivariate grid $\{(y_{j}, w_{q})\}_{j\in Z, q \in \mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}}\newcommand{\E}{\mathbb{E}}$ and at different points in time $t_{l+1}$ and $t_{l}$. It also states the rate of convergence of the algorithm. Define \mbox{$\Gamma_{1}(y_{i}, w_{p},t_{l}; y_{j}, t_{l+1})=\displaystyle\sum_{q=1}^{\infty}\Gamma_2(y_{i}, w_{p},t_{l}; y_{j}, w_{q}, t_{l+1})\sqrt{\Delta w}$}.
\begin{proposition}\label{prop:convHes}
Let $H(y,T)$ be such that $|H(y,T)-H(y',T)|<C|y-y'|$ for a positive constant $C$, and for $|y-y'|<\Delta y$. Let $v_{ip}(t_{l})$ be defined for a set of dates $\{t_{l}\}_{l=1,\dots,L}$, with $t_{L}=T$, as follows:
\begin{alignat}{4}\label{eq:recursiveHes}
	v_{ip} (t_{l}) &= \max\Bigl\{ H(y_{i},t_{l}), \sum_{j=1}^{\infty} \Gamma_{1}(y_{i}, w_{p},t_{l}; y_{j}, t_{l+1}) H(y_{j}, t_{l+1}) \sqrt{\Delta y} \Big\}, &	\ & \mathrm{for}	 & \ & l=L-1 ,	\\
	\label{eq:recursiveHes2}
	v_{ip} (t_{l}) &= \max\Big\{H(y_{i},t_{l}), \sum_{j,q=1}^{\infty}\mspace{-3mu}\Gamma_2(y_{i}, w_{p},t_{l}; y_{j}, w_{q}, t_{l+1})\mspace{1mu}v_{jq}(t_{l+1})\mspace{1mu}\sqrt{\Delta y \Delta w} \Big\}, & & \mathrm{for } & & l=1,\dots, L-2.
\end{alignat}

Then, for each $t_{l}$ in $\{t_{1},\dots, t_{L-1}\},$ the approximated values $v_{ip}(t_{l})$ defined in \eqref{eq:recursiveHes} and \eqref{eq:recursiveHes2} converge to the true value $V(y_{i}, w_{p}, t_{l})$ with an approximation error of the order $O\bigl(\underline\Delta^{2}\bigr)$, with $\underline \Delta = \sqrt{(\Delta y)^{2} + (\Delta w)^{2}}$.
\end{proposition}
\begin{proof}
See Section B of the supplementary online Appendix.
\end{proof}

For $t=t_{L-1}$,  Equation \eqref{eq:recursiveHes} gives the price of a European option in the Heston model.
Because the payoff function $H(y,T)$ only depends on the value $y$ taken by the underlying asset at $t_{L}=T$, the computed price $v_{ip}(t_{L-1})$ depends on the stochastic variance only through the conditioning value $\sigma^{2}_{t_{L-1}}=w_{p}$.
For this reason, we can use $\{ \Gamma_{1}(y_{i}, w_{p},t; y_{j}, T)\}_{j \in \mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}}\newcommand{\E}{\mathbb{E}}$ instead of $\{\Gamma_2(y_{i}, w_{p},t; y_{j}, w_{q}, T)\}_{j, q \in \mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}}\newcommand{\E}{\mathbb{E}}$.
Using  the values $\{ \Gamma_{1}(y_{i}, w_{p},t; y_{j}, T)\}_{j \in \mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\R}{\mathbb{R}}\newcommand{\E}{\mathbb{E}}$ in \eqref{eq:recursiveHes} is equivalent to applying the projection step on the Fourier transform $\hat{G}_{1}(y_{i}, w_{p}, t; \lambda,T) = \int \mspace{-3mu}d\kappa\mspace{1mu} \hat{G}_{2}(y_{i}, w_{p}, t; \lambda,\kappa, T)$. It is the univariate function $\hat{G}_{1}$ and not the bivariate  $\hat{G}_{2}$ that appears, for instance, in the original work by Heston (1993) for a European option.
Figure \ref{fig:projection-2d}  graphically presents the projection and recursion steps in the bivariate case.

	\medskip
	{\centering
	[Figure  \ref{fig:projection-2d} about here]

	}

In the implementation, we truncate the summations in \eqref{eq:recursiveHes} and \eqref{eq:recursiveHes2}, so that the grid $\{(y_{j},w_{q})\}_{j=1,\dots,N;q=1\dots, W}$ has $N\times W$ points.
The $N\times W$ matrix of computed prices at time $t=t_{l}$  is denoted by $\bm v_{2}(t_{l})$ , that is $\bm{v}_{2, jq}(t_{l})=v_{jq}(t_{l})$. Let $\mathbf{\Gamma}_{2}(y_{i}, w_{p}, t_{l}; t_{l+1})$ be the $N\times W$ matrix of the approximated transition probabilities from the initial point $(y_{i}, w_{p})$ to the end points of the entire grid $\{(y_{j},w_{q})\}_{j=1,\dots,N;q=1\dots, W}$. As in Section \ref{sec:example_european}, we integrate the normalization parameter $\sqrt{\Delta y \Delta w}$ in the definition of the transition matrix. We then have that $\mathbf{\Gamma}_{2, jq}(y_{i}, w_{p}, t_{l}; t_{l+1})=\Gamma_2(y_{i}, w_{p},t_{l}; y_{j}, w_{q}, t_{l+1})\sqrt{\Delta y \Delta w}$.
Let $\phi_{j}=\{\hat{e}_{j}(-\lambda_{r})\}_{r=1,\dots, R}$ and $\varphi_{q}=\{\hat{\varepsilon}_{q}(-\kappa_{z})\}_{z=1,\dots, Z}$ be the values of the functions $\hat e_{j}(-\lambda)$ and $\hat \varepsilon_{q}(-\kappa)$ sampled at the grids $\{\lambda_{r}\}_{r=1,\dots, R}$ and $\{\kappa_{z}\}_{z=1,\dots, Z}$, respectively.
Furthermore, we define the $R \times N$ matrix $\bm{\phi} = (\phi_{1}, \dots, \phi_{N})$, the
$Z \times W$ matrix $\bm{\varphi} = (\varphi_{1}, \dots, \varphi_{W})$, and the $R \times Z$ matrix $\mathbf{\hat G}_{2}(y_{i}, w_{p}, t_{l}; t_{l+1})$ with entries $\mathbf{\hat G}_{2, rz}(y_{i}, w_{p}, t_{l}; t_{l+1})=\hat G_{2}(y_{i}, w_{p}, t_{l}; \lambda_{r}, \kappa_{z}, t_{l+1})$.
Then, we can write the coefficients of the projection step \eqref{eq:projectionHesG}  in matrix form as:
\begin{equation}\label{eq:projectionHesGMatr}
	\mathbf{\Gamma}_{2}(y_{i}, w_{p}, t_{l}; t_{l+1}) = \bm\phi'  \mathbf{\hat G}_{2}(y_{i}, w_{p}, t_{l}; t_{l+1}) \bm\varphi\sqrt{\Delta y \Delta w}.
\end{equation}
The recursive step \eqref{eq:recursiveHes2} becomes the following:
\begin{align}\notag
	v_{ip}(t_{l}) = \max \Bigl\{H(y_i, t_l), & \sum_{j=1}^{N}\sum_{q=1}^{W}  \Gamma_2(y_{i}, w_{p},t_{l}; y_{j}, w_{q}, t_{l+1})v_{jq}(t_{l+1})\sqrt{\Delta y\Delta w}\,\Bigr\}\\
	\label{eq:recursiveHes2Matr}
	=\mspace{9mu} &\max \bigl\{H(y_i, t_l),  \mathbf{\Gamma}_{2}(y_{i}, w_{p}, t_{l}; t_{l+1}) \bm : \bm v_{2}(t_{l+1})\, \bigr\},
\end{align}
where the symbol ``$\bm :$'' denotes the Frobenius, or entry-wise, product.

In Section D of the supplementary online Appendix, we show how to speed up the computation of the matrices $\mathbf{\Gamma}_{2}(y_{i}, w_{p}, t_{l}; t_{l+1})$ by taking advantage of the space translation invariance property of transition densities.
Our method contains the Fast Fourier Transform (FFT) as a special case. In the FFT, the univariate grids for  $\lambda$ and $\kappa$ are automatically set, which can sometimes cause an imprecise reconstruction of the $\mathbf{\Gamma}_{2}(y_{i}, w_{p}, t_{l}; t_{l+1})$ matrices.

\subsection{Numerical illustrations in the Heston model}\label{sec:NumericalHeston}
We investigate the performance of our method in a standard affine model such as the \cite{heston1993closed} model.  We study an American option, written on an asset $S_t$, which pays discrete dividends and that evolves according to the following stochastic volatility model:
\begin{equation}\label{eq:stochvol}\begin{split}
        dX_{t} &= \Bigl(r -\frac{1}{2}\sigma_{t}^{2}\Bigr)dt + \sqrt{\sigma_{t}^{2}}\cdot dW_{1t},\\
        d\sigma_{t}^{2} &= \beta\bigl(\sigma_{LT}^2 - \sigma_{t}^{2}\bigr)dt + \omega \sqrt{\sigma_{t}^{2}}\cdot dW_{2t}, \quad E(dW_{1t}\cdot dW_{2t})=\rho dt.
\end{split}\end{equation}
In Equation \eqref{eq:stochvol}, $X_{t}=\log S_{t}$, and $\sigma_{t}^{2}$ is the variance process. We work with $X_t$ to be able to implement the space invariance property of the transition matrices, as outlined in Section D of the supplementary online Appendix. We conduct two simulation studies.
In the first, the American call has a time to maturity of one year, and $3$ dividends worth $d=2$ are distributed at $t_{h}=0.25,0.5,0.75$.
In the second, the time to maturity remains one year, but a single large dividend $d=10$ is paid out after six months.
The process parameter values are the following: $r=0.05$,  $\sigma_{LT}=0.2$, $\beta=2$ and $\omega=0.2$. Moreover we choose the parameter $\rho$ to be equal to zero. We compute the price for an at-the-money option ($S_{0}=K=100$).
The benchmark method in this analysis is a finite-difference (hereafter \emph{FD}) numerical solution of the partial derivatives equation (\textit{PDE}) that describes the evolution
of the price process $V_t$ of the American call. We implement an alternating direction implicit (\emph{ADI}) variant of the finite-difference scheme.
For a recent discussion of schemes similar to  \emph{FD}, see, for instance, \cite{in2010adi}. This implementation is equivalent to a Crank-Nicolson scheme, which in standard problems converges at a rate $O\bigl((\Delta t)^{2}\bigr)$, where $\Delta t$ is the temporal discretization interval.
In both the \emph{FD} scheme and the recursive projections, the evolution of the option price $V_t$ is studied on a rectangular grid in the space $(X,\sigma^{2})$, with $X \in [\log(K) - 10\sigma_{LT}\sqrt{T}, \log(K) + 10\sigma_{LT}\sqrt{T}]$ and $\sigma^{2} \in [0, 0.3]$. In the \emph{FD} scheme, the parameter $m_s$ gives the number of equally spaced grid points in the $X$ direction, and $m_v$ gives the number of equally spaced grid points in the $\sigma^{2}$ direction, so that the grid points are $\{(X_i, \sigma^2_p)\}_{i=1,\dots,m_s; p=1,\dots,m_v}$.  The parameter $L_{T}$ gives the number of time steps used.
In the recursive projections, under a sampling scheme we define $\Delta y=2^{-J}a$, where $a$ is a positive constant that gives the step of the $\{y_{j}\}_{j=1,\dots,N}$ grid when $J=0$. Describing the convergence of the recursive projections in terms of the parameter $J$ emphasizes how the approximation error decreases each time the number of grid points is doubled. Similarly, $\Delta w = 2^{-J_{w}}a_{w}$, where $a_{w}$ is the step with $J_{w}=0$ of the $\{w\}_{p=1,..,W}$ grid in which the $\sigma^{2}_{t}$ variable takes values.

Assuming the contemporaneous correlation $\rho=0$ simplifies the implementation of the \emph{FD} scheme, in the sense that neglecting the correlation between $X_{t}$ and $\sigma^{2}_{t}$ makes the \emph{FD} scheme easier to code and faster.
On the other hand, the speed and complexity of the recursive projection method are unaffected by the value chosen for the parameter $\rho$. The correlation is addressed in the Green function  $G_2(x, \sigma^{2}_{t}, t; y, w, T )$ and consequently in the coefficients of the matrix $\hat{\mathbf{G}}_{2}$. Because the speed of the method depends on the number of entries in the $\hat{\mathbf{G}}_{2}$ matrix, and not on the values taken by the entries, it is clear that the choice of $\rho$
does not affect the convergence rate of the recursive projections. This feature is the first advantage of the recursive projection over finite-difference schemes. This simulation study will then give a lower bound to the difference in speed between the recursive projections and the \emph{FD} scheme.
To price an American option on dividend-paying stocks, we should implement the \emph{FD} scheme-equivalent of the recombining tree. Doing so is practically unfeasible because it would mean computing at each ex-dividend date a new option price at each point of the grid.
Instead,  at each ex-dividend date $t_{h}$ and at each grid point $(X_{i}, \sigma^2_{p})$, we opt for comparing  the intrinsic value $H(X_{i}, t_h)$ with the continuation value $V(\widetilde{X}^{d}_{i}, \sigma^2_{p}, t_h)$, where $\widetilde{X}^{d}_{i}$ is the value of the $X$ grid closest to $\log(e^{X_{i}}-d)$.
This choice amounts to perturbating the \emph{FD} scheme at each ex-dividend rate, which could translate into a convergence slower than the theoretical $O\bigl((\Delta t)^{2}\bigr)$. This feature is a second advantage of the recursive projection over the finite-difference schemes, because, as we explained in Section \ref{sec:BSillustration}, the recursive projections can easily adapt to discrete dividends without their affecting the convergence properties of the algorithm.
The recursive projections achieve convergence quickly in the $\sigma^{2}$ direction. The method does not seem to improve by setting a resolution level greater than $J_{w}=4$; thus, we keep this value fixed throughout our simulations.
The \emph{FD} scheme is also not very sensitive to the number of points used in the $\sigma^{2}$ direction. We find no improvement beyond $m_{v}=31$.

Figure \ref{fig:Hes3divConv} shows the results for the 3-dividend case. The true value used to compute the pricing errors is $7.397$, obtained with the resolution level $J=13$. The graph on the right displays the pricing error of the \emph{FD} scheme as a function of the time discretization parameter $L_{T}$. Each line is relative to a different value of the spatial discretization parameter $m_{s}$. The time labels are all relative to the $m_{s}=3200$ curve.
The \emph{FD} scheme with $L_{T} =2048$ and $m_{s}=6400$ delivers a value within 1$bp$; thus, we assume that the methods have converged when the absolute value of the relative error is within 1$bp$ of $7.397$. The graph on the left plots the relative pricing error of the recursive projections against the resolution level $J$. The regression line on the left graph shows that the estimated slope is almost exactly the slope of $-2$ predicted by the theoretical convergence results of Proposition \ref{prop:convHes}.
 The \emph{FD} is at least one order of magnitude slower. Compare, for instance, the computation time needed to deliver a 4$bp$ error (2$s$ against 65$s$), or a 1$bp$ error (8$s$  against 130$s$). Figure \ref{fig:Hes1divConv} compares the convergence speed of the two methods in the 1-dividend case. The true value of 7.302 is obtained by the recursive projection method with $J=13$. The \emph{FD} scheme requires 48 seconds to reach a $5 bp$ relative error,  with parameters $m_{s}=400$ and $L_{T}=2048$. The bottom curve, relative to $m_{s}=200$, shows that the method does not converge for smaller values of the space discretization parameter.
The small $5 bp$ bias of the \emph{FD} is due to the large value of the dividend $d$ and the perturbation of the scheme at each dividend date. The rate as a function of the resolution level $J$ at which the recursive projections attain the $1 bp$ error band is approximately -2, as theoretically predicted.

The reason for the difference in speed between the recursive projections and the \emph{FD} scheme lies in the fundamentally different way  finite differences and quadrature methods deal with time stepping.
Both methods achieve time stepping through matrix multiplications. But while the number of time steps in the \emph{FD} is of the order of $2^{9}$ or higher, the recursive projections only need 3 or 4 time steps, one per divided payment, plus the expiry date.
The size of the parameters $L_{T}$, $m_{s}$ and $m_{v}$ determines the efficiency of the implementation of the \emph{FD} scheme.
If we compare the magnitude of the parameters $L_{T}$, $m_{s}$ and $m_{v}$ that we need to obtain convergence with the values of the equivalent parameters in \cite{in2010adi}, we find that our implementation is close to the most recent ones in the literature.
While specific implementations could marginally improve on ours, we think that we give a fair representation of the potential of the two techniques.
We remind that the computational time per time step is underestimated in our simulation, since the assumption of $\rho = 0$ reduces the number of intermediate steps in the $ADI$ implementation of the \emph{FD} scheme. Finally, if we include jumps in the process of the underlying stock, as we do in our empirical application, the numerical complexity of the recursive projections remains exactly the same as in the stochastic volatility case.
Introducing jumps in the underlying process while keeping the finite differences viable from a computational point of view asks for technical devices (see for instance \cite{d2005robust}) which are model specific and not yet implemented in conjunction with stochastic volatility\footnote{In Appendix H of the supplemental online Appendix we elaborate further on the implications of the choice $\rho=0$, on the relative advantages of recursive projections and $ADI$, and compare in detail our implementation with the one of \cite{d2005robust}.
}.

Another notable difference between the \emph{FD} and the recursive projection method is that the latter demands far fewer changes to adapt to different pricing problems.
In Equation \eqref{eq:matrix}, the matrix $\mathbf{G}(t;T)$ depends only on the dynamics of the underlying asset and not on the payoff. We can compute it once for all and use it to price different options with different payoffs, because the payoff functional form only impacts the vector $\mathbf{H}(T)$.  Such a design is particularly suited for object-oriented programming, which is often used in quant desks.
In finite-difference schemes, we cannot price options with different payoffs through the use of the same transition matrices, as boundary conditions affect the way the matrices are computed.

\section{Numerical applications and empirics}\label{sec:boundaries_empirics}
\subsection{Numerical comparison of early exercise boundaries}\label{sec:early_boundaries}

In this section, we compare the early exercise boundary implied by the Black-Scholes model with those implied by the Heston stochastic volatility and the Merton jump-diffusion models. We study two cases in which i) the stock distributes a continuous dividend yield and ii) the stock distributes a discrete dividend. Combining cases i) and ii) with the different modelling assumptions for the underlying asset and different maturities leads to very different patterns. For instance, in the discrete dividend case, the early exercise boundary is lower under the Black-Scholes model than under the Heston model, whereas in the continuous dividend case, the opposite is true. Hence, by modelling a discrete dividend as a continuous yield, we can draw misleading conclusions in an empirical evaluation of suboptimal non-exercise.
The exercise boundary $S^*_{t}$ for an American call with a continuous dividend yield is defined as the lowest value of $S_{t}$ such that $S_{t}-K \geq C(S_{t},T,K)$. If the value of the current stock is above $S^*_{t}$, then it is optimal for the call holder to exercise his option. With a discrete dividend, it can only be optimal to exercise the call option on the days immediately before the ex-dividend dates $t_{h}$. The exercise boundary $S^*_{t}$ for an American call option is then defined as the lowest value of $S_{t_{h}}$ such that $S_{t_{h}}-K \geq C(S_{t_{h}}-d,T,K)$.

	\medskip
	{\centering
	[Figure  \ref{Heston_exercise_boundary} about here]

	}

In Panel A of Figure \ref{Heston_exercise_boundary}, we plot the early exercise boundary for the Heston and Black-Scholes models for an American call option with a continuous dividend yield $r_d=0.03$ (right graph) and with an equivalent quarterly discrete dividend $d=1.38$ (left graph). We choose $d=1.38$ to have an equivalent total annual dividend between the continuous dividend yield $r_d=0.03$ and the discrete dividend case. Indeed, $1.38=0.03 S^*/4$, where $S^*=184$ is the critical stock price under the Black-Scholes model in the dividend yield case for maturity $T=0.5$. We use the following set of representative parameters: $T=0.5$, $K=100$, $r=0.05$, $\sigma_0=0.2$, $\omega=0.1$, $\sigma_{LT}=0.3$, $\beta=4$, and $\rho=-0.5$ (\cite{adolfsson2013representation}).
For comparison, we follow \cite{heston1993closed}, and we use the Black-Scholes model with a volatility parameter that matches the (square root of the) variance of the spot return over the life of the option in the Heston model.
When the stock distributes a regular quarterly dividend, there are only two dividend payments during the life time of an option with maturity $T=0.5$, and it is immediately before the payment dates that it can be optimal to exercise the option.
In our example, the two dividend payments occur immediately, at $t=0$, and at $t=0.25$,  corresponding to a time to maturity of $0.5$ and $0.25$, respectively.
At both dates, the value of the exercise boundary is lower under the Black-Scholes model. With a continuous dividend yield, the Heston early exercise boundary is always below the Black-Scholes boundary, whereas with discrete dividends, the opposite is the case.

Although the findings in the continuous dividend case are in line with those of \cite{adolfsson2013representation}, the findings in the discrete dividend case are entirely new.
This difference warrants further intuitive discussion.
Assume there is only one discrete dividend to be paid. The continuation value of the call option immediately after the ex-dividend date is that of a European call with the remaining time to maturity. When $\rho \leq 0$, the price of European options for a deep in-the-money call, where early exercise could be optimal, is higher in the Heston case than in the Black-Scholes case (see \cite{heston1993closed,hull1987pricing}).
For instance, in the left graph of Panel B of Figure \ref{Heston_exercise_boundary}, for a time to maturity of $0.25$, this would be the case in the range of stock prices of approximately $150$. Even by taking into account the dividend drop in computing the continuation value, the ex-dividend stock price should remain in the region where the Heston price is higher.
We can repeat the same argument for a number of discrete dividends sufficiently small (typically of the order of a couple percent) to prevent the stock price from falling in the price range where the call has more value under the Black-Scholes model. The behavior of the boundary with a continuous dividend is less straightforward to grasp.
Following \cite{Kim:1990aa} and \cite{jamshidian1992analysis}, we can decompose the value $V(S_{t},t)$ of an American option into two components, namely, the European value $V^{E}(S_{t},t)$ and the early exercise premium $V^{A}(S_{t},t)$, such that:
\begin{align}\label{eq:EarlyExPremium}
	V(S_{t},t) & = V^{E}(S_{t},t) + V^{A}(S_{t},t)	\\ \notag
			& =e^{-r(T-t)} \mathbb{E}\bigl[ (S_{T}-K)_{+}\big | S_{t}, \sigma^{2}_{t} \bigr] +\int_{t}^{T}\mspace{-9mu}e^{-r(s-t)} \mathbb{E}\bigl[ (r_{d}S_{s} -rK)\mathbb{I}_{(S_{s}>S^{*}_{s})}\big | S_{t}, \sigma^{2}_{t} \bigr]ds,
\end{align}
where $S^{*}_{s}$ is the early exercise boundary at time $s$ and $\mathbb{I}_{(S_{s}>S^{*}_{s})}$ equals one if, at time $s$, the stock is in the exercise region, otherwise zero. We can interpret $V^{A}(S_{t},t)$ as a  continuum of European call options with maturity $T-s$, strike price $S^{*}_{s}$, and payoff $r_{d}S_{s} -rK$.
For each of these European options, we can apply the results of Table II in \cite{hull1987pricing} who compare the values of
European options under general stochastic volatility dynamics with the Black-Scholes price.
Call values under the stochastic volatility assumption are lower when the contracts are at-the-money and $\rho \leq 0$.
The continuum of contracts composing the $V^{A}(S_{t},t)$ are at-the-money when $S_{s}=S^{*}_{s}$. As confirmed from our numerical simulations, the $S^{*}_{s}$ values are distributed in the region immediately above $S=150$, that is, exactly where the price of the American option under the Heston model is lower than that under the Black-Scholes model, and explain the negative bump in the right graph of Panel B of Figure \ref{Heston_exercise_boundary}.

Similarly, we can characterize the early exercise under the Merton jump-diffusion model, where the asset $S_{t}$ evolves according to the following jump-diffusion process:
\begin{equation}
    \frac{dS_t}{S_t}=(r-r_d-\gamma \upsilon)dt + \sigma_M dW_t +(\psi-1) dq_t,\label{eq:Merton}
\end{equation}
 where $r_d$ is the continuous dividend yield paid by the asset, and $\sigma_M^2$ is the instantaneous variance of the return conditional on no jump arrivals. The Poisson process, $q(t)$, is independent of $W_{t}$, and such that there is a probability $\gamma dt$ that a jump occurs in $dt$, and $1-\gamma$ probability that no jump occurs.
The parameter $\gamma$ represents the mean number of jumps per unit of time. The random variable $\psi$ is such that $\psi-1$ describes the percentage change in the stock price if the Poisson event occurs, and $\upsilon=\mathbb{E}[\psi-1]$ is the mean jump size.
We further make the standard assumption (for instance, see \cite{Amin:1993a,bakshi1997empirical}) that $\log(\psi) \sim N(\mu_{\psi},\sigma_{\psi}^2)$. If $\gamma=0$, then we recover the standard Black-Scholes model with no jumps.
We use the following set of representative parameters: $K=40$, $r=0.08$, $\gamma=5$, $\sigma_M^2=0.05$, $\sigma_{\psi}^2=0.05$, $\mu_{\psi}=0$ (\cite{Amin:1993a}).
We set the volatility parameter in the Black-Scholes model  equal to the volatility of the underlying return over the life of the option in the Merton model.

	\medskip
	{\centering
	[Figure  \ref{Merton_exercise_boundary} about here]

	}

In the Panel A of Figure \ref{Merton_exercise_boundary}, we plot the early exercise boundary for the Merton and Black-Scholes models for an American call option with a continuous dividend yield $r_d=0.05$ (right graph) and in the case in which the stock pays an equivalent quarterly discrete dividend \mbox{$d=1.125 $}~\footnote{As before, we take $d=1.125$ because $1.125=0.05 S^*/4$, where $S^*=90$ is the critical stock price with the dividend yield $r_d=0.05$ for maturity $T=0.5$.} (left graph). As for the Heston case, the results on the continuous dividend case are in line with the existing literature, e.~g., \cite{Amin:1993a}, and we provide new insights into the discrete dividend case.
To interpret the graphs in Figure \ref{Merton_exercise_boundary}, we have to make an important distinction. For short maturity options, the jump component in Equation \eqref{eq:Merton} dominates the diffusion component. As explained in \cite{Amin:1993a} and \cite{merton1976option}, the result is higher prices for short maturities under the Merton model than under the Black-Scholes model. We call this effect the jump effect. For longer maturities, the jump effect no longer dominates the diffusion component but instead creates an interplay that makes the jump-diffusion process observationally similar to a stochastic volatility process. For a discrete dividend, both the jump effect and the stochastic volatility effect, as previously discussed in the Heston case, predict a higher boundary in the Merton case than in the Black-Scholes case, which holds true for all maturities. This result is exactly what we find in the left graph of Panel A of Figure \ref{Merton_exercise_boundary}.
For a continuous dividend case, we have that the jump effect dominates for short maturities, giving a higher boundary under the Merton model than under the Black-Scholes model. For longer maturities, the jump effect diminishes, and the boundary behaves as in the stochastic volatility model, that is, taking lower values than in the Black-Scholes case. These insights explain the crossing of the early exercise boundary that we observe in the right graph of Panel A, Figure \ref{Merton_exercise_boundary}.

A key numerical finding of this section is that the early exercise is more likely under the Black-Scholes model when we face discrete dividends. In the next section, we assess the empirical consequences of this finding for the cost of suboptimal non-exercise.

\subsection{Empirical analysis}\label{sec:empirics}
In this section, we apply the recursive projection method to characterize the early exercise boundary of a large sample of call options. The options have a maturity of less than six months and are written on dividend-paying stocks, which are part of the Dow Jones Industrial Average Index (DJIA). The sample comprises daily observations between January 1996 and December 2012. We investigate the early exercise decision of call holders in light of the different values
that the exercise boundary can take under distinct modelling assumptions for the underlying asset.  Following the procedure suggested by \cite{pool2008failure}, we first check which contracts should be exercised by comparing the intrinsic value immediately before the dividend payment with the continuation value on the ex-dividend day. We quantify how much is economically lost in the case of a suboptimal non-exercise decision.
This amount depends on the continuation value and is model-specific. We compare the results obtained under three modelling environments, namely,
Black-Scholes, Merton jump-diffusion, and Merton jump-diffusion with the stochastic volatility dynamics of the Heston model.
\cite{Bates1996} was the first to suggest combining the Merton and Heston models, and therefore, we refer to this process specification as the Bates model.
Finally, whenever we find evidence of a suboptimal non-exercise decision, we show that trading costs alone cannot justify the behavior of investors.

In our empirical analysis, we price by fully taking into account the discrete nature of the dividend distributed by the underlying stocks and the American style of the call options, and we do so for the three pricing models. This feature is a peculiarity of our work, given that the standard empirical literature on options mainly focuses on European $S\&P 500$ options with a dividend yield or limits itself to the Black-Scholes model. We have already shown in an example in Section \ref{sec:BSillustration} how neglecting the discrete cash flow feature of the dividend payment, and approximating it with a continuous dividend yield, leads to a pricing error of $169bp$. In our empirical analysis, we also need to correctly take into account  the time of payment. A popular approach, when dealing with American options on dividend-paying stocks, is the so called ``escrow dividend'' model. Under this approximation, the option is priced as if it were a European contract, valued at the prevailing stock price minus the present value of all the dividends to be distributed during its remaining life. This technique correctly models that a long holder of a call option is unprotected from dividend distributions, but does not properly integrate the early exercise premium of an American option, and tends to underprice the option. As a consequence, an investor following this approach, could underestimate the value of the early exercise boundary, and exercise a contract that he would be better off by keeping alive. For example, on  May 10th, 2006, the Dupont stock closes at $45.71$ dollars. The call option with $K=30$ and $T=0.45$ should not be exercised if the continuation value is computed correctly, but the option should be exercised if the continuation value is approximated with a European price. In this case, if an investor exercises his option wrongly, he suffers a loss of $200 bp$ if the spot price is modelled under the Merton dynamics\footnote{We elaborate further on this example in Section G of the supplementary online Appendix where we show that a correct modelling of the dividend as a discrete cash flow does matter when measuring the cost of sub-optimality in the empirics of the next section.}.

 Our choice for alternative modelling environments follows the empirical findings of \cite{bakshi1997empirical}, who suggest that jumps and stochastic volatility play a dominant role in pricing short-term options whereas modelling stochastic interest rates does not seem to significantly improve the pricing performance.
In addition, the choice of the jump arrival distribution in Equation \eqref{eq:Merton} is motivated by the work of \cite{bajgrowicz2016jumps}, who show that high-frequency data on individual stocks support the hypothesis that jump arrivals follow a simple low-intensity Poisson process.
Their findings also support the assumption in Merton (1976) that the jump component is nonsystematic, i.e., diversifiable, because of the absence of co-jumps affecting all stocks. This explains our choice of the Bates model, where $X_{t}=\log(S_{t})$, with:
\begin{align}\label{eq:MRT+HST}
    & dX_t =(r-r_d-\gamma \upsilon-\frac{1}{2}\sigma^2_{t})dt + \sigma_{t} dW_{1,t} +log(\psi) dq_t,\\ \nonumber
    & d\sigma_{t}^{2} =\beta \bigl(\sigma_{LT}^{2} - \sigma_{ t}^{2}\bigr)dt + \omega \sqrt{\sigma_{t}^{2}}\cdot dW_{2,t}, \quad E(dW_{1,t}\cdot dW_{2,t})=\rho dt,
\end{align}
and $\log(\psi) \sim N(\mu_{\psi},\sigma^2_{\psi})$.
To implement the recursive projections in the Bates model, we need the Green function implied by Equation \eqref{eq:MRT+HST}. Because the jump process is independent of the Brownian motions $ dW_{1,t}$ and $ dW_{2,t}$, we use the property that the characteristic function of a bivariate process is the product of the characteristic functions of the independent univariate processes. The transition matrices are then obtained as in Section \ref{sec:RecursHes}.
The Bates model shows how straightforward it is to adapt the recursive projections to more complex models.
The number of entries in the transition matrices for the Bates model is exactly the same as it is as for the Heston model.
In the calibrations, we have compared the computation of the transition matrices with both the FFT and the full sampling method outlined in Section \ref{sec:ProjectionHes}. Because we have found that the FFT correctly reproduces the $\mathbf{\Gamma}_{2}(y_{i}, w_{p}, t_{l}; t_{l+1})$ matrices, we have opted for using the FFT throughout the empirical exercise to take advantage of the further increase in speed provided by the FFT algorithm.

The daily data on all option attributes, the stock price, and the dividend distribution details are from Optionmetrics. We obtain the daily data on the interest rates from the Treasury constant maturities of the H15 report of the Federal Reserve\footnote{We access the Optionmetrics and the H15 databases through the Wharton Research Data Services (WRDS) research platform.}.
A total of 101,295 series of short-term options written on 30 stocks enter our database. The total number of records is approximately 9.5 million.
This number stresses the importance of a fast and versatile numerical method.
Table \ref{Table:Db_description} reports the number of quotes for each stock with a breakdown by maturity and moneyness. Our study focuses on the early exercise behavior of investors; hence, we focus on the in-the-money options, which are the category of options for which the number of quotes is the highest.

	\medskip
	{\centering
	[Tables  \ref{Table:Db_description} and \ref{Table:Param} about here]

	}

\subsubsection{Estimating the cost of suboptimal non-exercise}\label{sec:cost}

Table \ref{Table:Param} reports the results of the calibration for the three modelling frameworks. We obtain the parameters through the minimization of the implied volatility mean squared error, as in \cite{christoffersen2004importance}\footnote{ Section F of the supplementary online Appendix gives a detailed description of the data and the calibration procedure.}.
The first line of Table  \ref{Table:Param} displays the average values of the parameters calibrated on our sample of single stocks, whereas the second line reports the average values that \cite{bakshi1997empirical} obtain for the same parameters from contracts written on the  $S\&P 500$ index.
The parameters that rule the level of the volatility smile, namely, the Black-Scholes volatility $\sigma_{BS}$, the long term volatility $\sigma_{LT}$, and the spot volatility $\sigma_0$, are much higher in our single stock calibration than in the index calibration, which reflects the well known fact that an index is less volatile than its components.
Indeed, in our sample, the average Black-Scholes volatility is $29\%$, $\sigma_0$ is $28\%$, and the average long-term implied volatility is $32\%$, whereas for the index options, the same parameters take the values of $18.15\%$, $20\%$, and again $20\%$.
The jump parameters in the Bates model show that jumps are on average less frequent in the single stock case than in the index case ($\gamma_{Stocks}=0.5$ against $\gamma_{SP500}=0.61$), but the amplitude and variability are higher for single stocks ($\mu_{\psi,Stocks}=-0.12$ and $\sigma_{\psi,Stocks}=0.18$ respectively) than for the index  ($\mu_{\psi,SP500}=-0.09$ and $\sigma_{\psi,SP500}=0.14$). Given that the index is a diversified portfolio, it displays a jump whenever one of its constituents jumps, but the non-systematic nature of the jumps of single stocks (see \cite{bajgrowicz2016jumps})  attenuates the variability and magnitude at an aggregate level.
The remaining two parameters of the stochastic volatility component of the Bates model,  the correlation parameter $\rho$ and the volatility of volatility $\omega$, have a specific impact on the shape of the implied volatility smile (\cite{hagan2002managing,west2005calibration}). A negative $\rho$ implies a negatively sloped smile. The correlation parameter is in absolute value lower in the single stock case ($\rho_{SP500}=-0.52$ against $\rho_{Stocks}=-0.35$), meaning that the
implied volatility smile for the index is more negatively sloped than for individual stocks.
This finding is consistent with the findings of  \cite{bakshi2003stock} who describe the same relationship between the slopes of the index and of the individual stocks implied volatility smiles. \cite{bollen2004does} also find the same pattern, and explain it by relating the slope of the index smile to the buying pressure for index puts, with the demand for call options driving the shape of the smile of single stocks. The volatility of volatility $\omega$ determines the convexity of the implied volatility smile. The difference in the values taken by $\omega$ is striking. We find $\omega$ to be $75\%$ for stock options, whereas \cite{bakshi1997empirical} find a much smaller value of $40\%$ for short-term index options. This difference is due to the higher convexity of the implied volatility smiles of stock options versus that of index options, another feature also documented in \cite{bollen2004does}. The parameters $\rho$ and $\omega$  are related to the smile shape through the higher moments of the distributions of the returns of the underlying. A more negative $\rho$ generates a more negatively skewed distribution of index returns with respect to stock returns, whereas a higher $\omega $ in the single stock returns leads to a higher kurtosis than in the index return distribution.

In Table I of the supplementary online Appendix, we provide more detailed results, including a breakdown of the calibration by stock, and we show that the values of the calibrated parameters are homogeneous across stocks.

	\medskip
	{\centering
	[Table  \ref{Table:Param} about here]

	}

After having calibrated the models, we are able to compute the price $C(S_{t_{h}-1}-d, K, T)$ on the day previous to the dividend payment date $t_{h}$, by using the recursive projections. The price $C(S_{t_{h}-1}-d, K, T)$ is the continuation value of the option at date $t_{h}$, when the dividend $d$ is distributed. By comparing it with the intrinsic value $S_{t_{h}-1}-K$, we can assess which options should be exercised on $t_{h}-1$. If an option should be exercised (i.e., $C(S_{t_{h}-1}-d,K,T) \leq S_{t_{h}-1}-K$), then a positive open interest at the end of the day before ex-dividend ($OI_{t_{h}-1}>0$) measures the failure of investors to exercise the option. In this case, we calculate the suboptimal non-exercise percentage as the following ratio:
\begin{equation}\label{eq:percentage_unexercise}
    NE_{\%}=\frac{OI_{t_{h}-1}}{OI_{t_{h}-2}},
\end{equation}
i.e., the number of contracts outstanding at the end of the day $t_{h}-1$ to the total number of contracts outstanding at the end of day $t_{h}-2$.
The quantity defined in Equation \eqref{eq:percentage_unexercise} is an approximation of the actual non-exercise ratio, because it neglects a possible issue of new contracts on date $t_{h}-1$. This event is unlikely; indeed \cite{pool2008failure} test the approximation on a subsample of contracts for which they have the real exercise data. They conclude that the approximation (\ref{eq:percentage_unexercise}) is a precise description of the actual exercise behavior of option investors. The total amount of money that is left on the table due to suboptimal
non-exercise is given by the following formula:
\begin{equation}\label{eq:TotalLoss}
TL=100 \times OI_{t-1} \times [(S_{t-1}-K)-C(S_{t-1}-d,K,T)].
\end{equation}
The continuation value $C(S_{t-1}-d,K,T)$ depends on the model used for pricing; hence, the total loss due to suboptimal non-exercise (TL) is itself model-specific.

Table \ref{table:NII_lost} presents the results on the suboptimal non-exercise behavior of investors.

	\medskip
	{\centering
	[Table  \ref{table:NII_lost} about here]

	}

\noindent Table \ref{table:NII_lost} clearly shows that the optimal early exercise decision depends on the model used for the stock price. Under the Black-Scholes model, approximately $9.5\%$ of the outstanding contracts should be exercised, and the percentage decreases (approximately $7.5\%$) under the alternative models.
This result is consistent with the numerical findings of Section \ref{sec:early_boundaries},
where we show that, in the case of discrete dividends, the early exercise boundary under the Black-Scholes model is lower compared to that implied by the Merton and Heston models. As a general rule, contracts that should be exercised under the Merton or Bates models should also be exercised under the Black-Scholes model. In our sample, we find some exceptions to this rule because,
in Section \ref{sec:early_boundaries},
we choose the model parameters such that the total variance of the returns over the life of the option is the same in all models, whereas in real data, this condition may not hold. To give some examples, $4680$ options should be exercised under Black-Scholes but not under the Bates model, whereas the reverse is true only for $249$ contracts.
Similarly, we find that $2872$ options should be exercised under the Black-Scholes model but not under the Merton model, whereas the opposite occurs with only $53$ options. The first important lesson we learn is that, by allowing for more sophisticated models than the Black-Scholes model, the number of contracts that should be optimally exercised decreases by almost $25\%$. The suboptimality figures are model-dependent and may be a consequence of the calibration procedure.
The comparison between our calibration results and the results of \cite{bakshi1997empirical} are reassuring in terms of the reliability of our calibration method. To justify the suboptimal behavior found in our sample, we should obtain unreasonably high values for the jump and intensity parameters.

A second piece of evidence that stands out from Table \ref{table:NII_lost} is that the percentage of investors who leave the options suboptimally non-exercised is higher under the Black-Scholes model than under the other models, $39\%$ versus approximately $30\%$. We compute these percentages in accordance with Definition (\ref{eq:percentage_unexercise}).
If we restrict our attention to the $1965$ contracts in our sample that should be exercised under the Black-Scholes model but not under the Merton model or the Bates model, we find a striking $81\%$ of no-exercise. These results may suggest that investors do not limit themselves to a Black-Scholes world when evaluating their options but rely on more sophisticated models that include jumps or stochastic volatility.
Even if this evidence is a considerable step towards understanding the investor decision-making process, it does not fully solve the puzzle. Indeed, even in the Merton and Bates models, we still find a high percentage of suboptimal non-exercises, which leads to a global loss of approximately $130-140$ million dollars, down approximately
$30\%$ from the loss of $206$ million dollars in the Black-Scholes model.

A second possible explanation of the early exercise puzzle is that investors wait until the options are deeper in-the-money before exercising. If we restrict ourselves to the options that should be exercised only under the Black-Scholes model but not under the Merton or Bates models, then the average moneyness of the contracts, measured by the delta of the options, is $0.962$. If we consider all the contracts that should be exercised under the Black-Scholes model, including those that should also be exercised under the other two model specifications, the average moneyness rises to $0.986$.
Figure \ref{Figure:moneyness} shows how relevant the moneyness of the contract is in the non-exercise decision of investors: the more in-the-money the option is, the smaller the number of suboptimal non-exercises.
Investors may not respond immediately to favorable stock price movements and may take some time before reacting and optimally exercising their option, which would be in line with the behavior on the early prepayment of mortgages documented by \cite{stanton1995rational}.

	{\centering
	[Figure \ref{Figure:moneyness} about here]

	}

Our summary results regarding the exercise decision in the Black-Scholes model are in line with those obtained by \cite{pool2008failure}.
In that work, the authors apply the early exercise decision rule to all options series by using the Black-Scholes model with historical volatility and find that $53.1\%$ of investors leave their options unexercised when instead they should have been exercised. Their data span over ten years (from 1996 to 2006) and to compare our results with theirs, we divide our sample into two subsamples, the first spanning the years 1996-2006 and the second spanning the years 2006-2012. Then, we calculate the average percentage of suboptimal non-exercise in the two subsamples and find that the percentage of suboptimal non-exercise under the Black-Scholes model is approximately $47\%$ in the first subsample and $37\%$ in the second. The decrease in the non-exercise behavior with time intimates that investors become more attentive in monitoring their investments. There is a small difference between our results ($47\%$) and the $53.1\%$ found in
\cite{pool2008failure}.  The explanation is most likely our focus on the constituents of the Dow Jones Industrial Average, whereas \cite{pool2008failure} consider all option series. It is likely that, for large-cap companies, stock and option prices are monitored more closely than they are on average.

Throughout our empirical investigation, we choose a model-based approach to calculate the continuation value of the option $C(S_{t_{h}-1}-d,K,T)$. We could have also used a market-based approach where the continuation value is the market price of the option. The market-based approach checks whether the quantity $C_{MKT}(S_{t_{h}-1}, K, T)-(S_{t_{h}-1}-K)_+$ equals $0$, where $C_{MKT}(S_{t_{h}-1}, K, T)$ is the observed market price at $t=t_h-1$. As discussed in \cite{pool2008failure} and in \cite{barraclough2012early}, the market-based approach has shortcomings. The most important is that it does not make it possible to calculate the total loss due to suboptimal non-exercise, which we do in Equation (\ref{eq:TotalLoss}). In addition, the bid-ask spread and the discreteness of the prices make it difficult to decide which $C_{MKT}$ should be used. For all of these reasons, we follow \cite{pool2008failure} and \cite{barraclough2012early} and use a model-based approach with Equation \eqref{eq:percentage_unexercise} to account for the actual exercise behavior of options investors. \cite{barraclough2012early} only use the market-based approach as a useful model-free test to confirm the presence of suboptimal non-exercise behavior. They find that the market-based approach gives a magnitude of suboptimality that is comparable to that implied by the model-based approach. This last piece of evidence is an additional argument against the possible objection that an incorrect calibration of the model parameters may be the source of the suboptimal exercise figures.

\subsubsection{The role of fees}

According to the recent literature on option prices (\cite{jensen2016early, christoffersen2014illiquidity}), trading costs and financial frictions in general strongly affect both the option prices and the early exercise decision of American options. In this section, we investigate whether the suboptimal non-exercise behavior of investors is due to the trading costs that investors face when exercising their options.

Following \cite{pool2008failure}, we model the costs of exiting a long call position as a per share lump sum $\mathcal{F}$ that the investor must pay at the moment he decides to exercise. The specific value of $\mathcal{F}$ depends on how the exit is accomplished according to the different possible objectives of the investor. The most expensive value of the fee $\mathcal{F}$ is attained when the investor wants to exercise the option and reenter into the same call position. \cite{pool2008failure} estimate an average value for the rollover costs $\mathcal{F}$ by using the commissions of the high-cost brokers, and they obtain a very conservative amount of $\mathcal{F}=0.4446$ dollar per share. A detailed description of the components of the fee $\mathcal{F}$ can be found in \cite{pool2008failure}.

To understand the role of the fees in the early exercise decision, we perform two different empirical exercises.
As a first check, we re-perform the exercise of Section \ref{sec:cost} and compute the loss due to a suboptimal non-exercise decision, but this time using $C(S_{t_{h}-1}-d,K+\mathcal{F},T) $ as the continuation value, and $(S_{t_{h}-1}-K-\mathcal{F}) $ as the intrinsic value. The fee value  is  $\mathcal{F}=0.4446$.
The fee $\mathcal{F}$ enters both in the exercise proceeds and in the continuation value. Indeed, at the moment of the exercise decision, the investor should decide whether to exercise and hence pay the exercise fee immediately or not exercise and postpone the payment of the exercise fee to a future date.
Accordingly, the calculation of the total amount of money that is lost due to suboptimal non-exercise is given by the following formula:
\begin{equation}\label{eq:TotalLossFee}
TL_{\mathcal{F}}=100 \times OI_{t-1} \times [(S_{t-1}-K-\mathcal{F})-C(S_{t-1}-d,K+\mathcal{F},T)].
\end{equation}
The second column in Table \ref{table:NII_lost} shows the summary results including the fee. They are not very different from those obtained without considering the fee (first column of Table \ref{table:NII_lost}). We can conclude that the inclusion of trading costs does not change the big picture on the suboptimal non-exercise of investors, as outlined in the previous paragraph.

As a second empirical exercise, we calculate the value of the fee that would justify the non-exercise decision of investors to detect possible additional costs that are not taken into consideration in the fee  $\mathcal{F}$. To do so, for each option for which $C(S_{t_{h}-1}-d,K,T)<(S_{t_{h}-1}-K)$, but that is not optimally exercised by some of the investors, we compute the value of the implied fee $\mathcal{F}$ that would justify the non-exercise decision. It amounts to numerically finding the zero of the following function:
\begin{equation}
    f(\mathcal{F})=C(S_{{t_{h}-1}}-d,K+\mathcal{F},T)-(S_{t_{h}-1}-K-\mathcal{F}).
\end{equation}
The results are reported in Table \ref{table:Fee}. The average implied fee is between $7$ and $8$ dollars per share, an incredibly high amount compared to the already conservative fee of $0.4446$ dollar per share estimated by \cite{pool2008failure}. No realistic hidden fees can sum up to 7 dollars per share, and the trading costs of exiting a long call option position cannot fully justify the suboptimal non-exercise behavior of investors.
We can interpret the difference between the implied fee of 8 dollars and the conservative fee of 0.4446 as an implied opportunity cost for the holder of the option to monitor the optimal exercise of the American option. The holder of the option chooses to spend this amount on alternative activities.

	{\centering
	[Figure \ref{table:Fee} about here]

	}

\section{Concluding remarks}\label{sec:conclusion}

We investigate the exercise behaviour of investors in a large database of 101,295 series of short-term American call options. \cite{pool2008failure} find that more than 40\% of the investors fail to optimally exercise their contracts.
We extend their analysis by including stochastic volatility and jumps to the process of the underlying stock.  In order to deal with the large option database and the repeated calculations required for the calibration and pricing, we develop an option pricing technique which is at the same time fast, precise and which can handle both multidimensional dynamics and cash dividend distributions. We start from the observation that, by monitoring the value of an option at discrete times and by sampling the value function of the contract on a finite grid of values of the underlying asset, we can describe the evolution of the price process in terms of elementary matrix operators.
The interpretation of the elements of such a matrix in terms of a functional projection allows us to extend the matrix approach to the pricing of contracts written on assets following processes whose transition probabilities do not have a known analytical expression in the direct space. The recursive projection method owes its speed to the simplicity of its algorithm, which is based on sampling. Moreover, our approach allows us to derive transition matrices that relate the option prices at points in time that are arbitrarily distant.
The number of time steps required is solely driven by features of the contracts, such as dates at which the contract needs to be monitored, and this is the main feature that makes our method faster than existing alternatives.

By applying our technique to the dataset, we can explain up to 25\% of the gain forgone due to suboptimal exercise decisions, as computed in \cite{pool2008failure}.
This result confirms the insights we obtain from the theoretical part of the paper. Indeed, we show that the exercise boundary is higher under the Merton and Heston models than under the Black-Scholes model if the dividend is discrete.
This result underlines the importance of the correct modelling of the dividend distribution.
We show that by modelling the dividend as a continuous yield instead of a discrete cash flow, the exercise frontier in the Bates model would have been lower instead of higher, and the suboptimal exercise behavior of  \cite{pool2008failure} would have been reinforced instead of mitigated.
We further try to check whether we can explain the remaining 75\% of suboptimal behavior in terms of transaction costs (\cite{jensen2016early}). We show that hidden transaction costs would need to be unrealistically large to explain the entire amount foregone by investors.
This observation leads us to interpret the implied transaction fee as a monitoring cost.

We can envisage two lines of research for further investigation: other possible applications of our method and improving the algorithm.
American option pricing is only a particular case of stochastic optimal control problems. We can think of applying the recursive projection method to other problems, such as the optimal portfolio allocation involving complex and path-dependent financial assets.
Currently, these types of complex problems are solved by using Monte Carlo simulations  (\cite{Detemple2003monte}) and our method could offer a more efficient computational alternative. Similarly our method makes possible other types of analysis with options based on individual stocks paying discrete dividends and sophisticated pricing models such as, e.g., the recent studies of tail risk (\cite{andersen2015risk, Andersen2016aa}), political uncertainty (\cite{Kelly2016a}), sector-wide government guarantees (\cite{kelly2016price}), and factor structure (\cite{christoffersen2015factor}). We can also use our method to extend the optimal control problem of \cite{hodder2007} to more complex dynamics for the decision variable.
The second possible extension of our work is the use of different function bases in the projection step. Candidate basis systems are functions that share the localization property with the indicator function but that may display higher regularity and enhance the convergence speed.
Faster convergence should not come at the cost of a more complicated numerical implementation, and therefore, research should go in the direction of functional projections that result in linear transformations of the sampled values (\cite{Sweldens1996lifting, Sweldens1998lifting}).

redpajama_set_name: RedPajamaArXiv