Numerical methods to value an option including risk aversion with a constant relative risk aversion function

Since the seminal work of Black and Scholes (1973), later extended by Merton (1973), the valuation of European options has relied on the assumption that the underlying asset follows a (GBM). Their framework provided the first closed-form solution and established risk-neutral pricing as a cornerstone of modern financial economics. Building on these theoretical foundations, various numerical estimation methods emerged, such as the finite-difference method, Adomian decomposition method and CRR proposed by Cox et al. (1979). Monte Carlo simulation techniques, such as the least squares method (LSM) proposed by Longstaff and Schwartz (2001), have been widely used to value American and Bermudan options with early-exercise features (Gapeev and Li, 2022). As Stentoft (2004) showed, LSM converges to the true expectation but may fail to converge to the real value of the option. Later, Létourneau and Stentoft (2014) modified the ordinary least square (OLS) regression by imposing constraints to mitigate the bias from the number of regressors, although asymptotic convergence was not established. More recently, Fabozzi et al. (2017) proposed weighted OLS estimations to address heteroscedasticity. Among these approaches, the CRR method remains the most widely used due to its simplicity, ease of construction and ability to replicate the discrete-time dynamics of the price of the underlying asset while providing tractable derivative valuation (Cox et al., 1979).

Discrete-time extensions have further advanced option valuation. Marín-Sánchez (2010) generalized the CRR through the MBRT method, while Milanesi (2020) incorporated isoelastic utilities and time-varying volatility in a trinomial grid for real options. More recently, Marín-Sánchez et al. (2021) employed quadrinomial trees to capture stochastic volatility, reinforcing the flexibility of lattice methods to manage nonlinear dynamics. These contributions highlight the adaptability of discrete models, yet they continue to assume risk-neutral investors.

The persistence of risk neutrality underscores a key gap in the literature: most frameworks do not explicitly account for heterogeneous investor risk preferences, despite empirical evidence that option prices embed valuable information about market risk preferences (Bliss and Panigirtzoglou, 2004). Existing attempts to incorporate preferences often rely on continuous-time models with high computational cost or lack direct adaptation to discrete-time settings.

Our contribution aims to fill this gap by proposing a discrete-time option pricing model that directly incorporates investor risk preferences through a CRRA utility specification. We adapt and extend the MBRT framework to value financial options using certainty-equivalent (CE) approach. Under the CE framework introduced by Beedles (1978), the value of a risky payoff is defined by its utility-adjusted equivalence to a certain outcome. Unlike classical models, our approach does not impose a full risk neutrality but instead reflects the actual risk attitudes of economic agents.

Importantly, while discounting remains risk-neutral, payoffs are adjusted using a CRRA utility function. This separation preserves flexibility and tractability: preferences need only satisfy standard monotonicity (more is preferred to less), while allowing for heterogeneous degrees of risk aversion. In the limiting case where γ = 0, our framework collapses to the risk-neutral benchmark, recovering classical pricing outcomes in both payoff equivalence and the discounting rule. In this sense, our framework generalizes traditional valuation models – including Black and Scholes (B-S), CRR and LSM – by preserving the risk-neutral discounting of expected payoffs (Björk, 2009), while introducing a nonlinear, preference-sensitive transformation of option values.

The remainder of the paper is organized as follows. Section 2 reviews the theoretical foundation of utility functions, focusing on CRRA and its link to consumption in option pricing. Section 3 derives the main equations, establishes the first two moments of the stochastic process, and details the MBRT framework with supporting algebraic expressions. Section 4 presents numerical experiments analyzing the impact of risk preferences on option values. Section 5 discusses the model's contribution to nonlinear risk-adjusted pricing, its practical relevance and limitations and suggests directions for future research. Section 6 concludes the study by emphasizing the importance of incorporating risk preferences into option valuation to enrich both theory and practice.

Decision-making under uncertainty is inherently challenging, requiring methodologies that extend beyond traditional approaches. In academia, risk is commonly modeled using expected utility theory (EUT), which reveals investor preferences by assigning values to the utility generated by alternatives choices, such as investment and consumption. Rooted in neoclassical finance (Shefrin, 2010), EUT follows a rigorous and systematic process that excludes behavioral assumptions. Kahnemann and Tversky (1979) argued that, from the perspective of EUT, decision analysis under risk constitutes a normative model of rational selection, applicable to a wide range of economic behaviors. A detailed review of utility functions and empirical methods for estimating risk aversion is provided by Chávez et al. (2017), who highlight key theoretical foundations and estimation techniques applied in financial decision-making contexts.

Risk aversion is generally quantified through two principal measures. The first is absolute risk aversion (ARA), also known as the Arrow–Pratt coefficient, expressed as $rA(x)=u″(x)/u′(x)$⁠, where $u(x)$ is the utility function, $u′(x)$ is the marginal utility of wealth and $u″(x)$ is its second derivative. To normalize by wealth, the relative risk aversion (RRA) coefficient is defined as $rR(x)=−x·u″(x)/u′(x)$⁠.

Consistent with both theory and empirical evidence, we adopt the CRRA utility function to model investor behavior under uncertainty. CRRA satisfies the following key properties: decreasing ARA, constant RRA and scale invariance. It provides analytical while capturing varying degrees of concavity in preferences. In the context of option pricing, CRRA facilitates both the adjustment of investor preferences and correct mathematical manipulation (Bliss and Panigirtzoglou, 2004; Brenner, 2015; Zhou, 2021). Experimental designs based on random lottery pairs further validate its empirical relevance, allowing estimation of individual risk-aversion parameters in controlled environments (Pareja-Vasseur and Baena, 2018).

The CRRA utility function takes the following form:

Here, $x>0$ represents wealth or a monetary payoff, and $γ∈R$ is the coefficient of CRRA utility function. This parameter defines attitudes toward risk:

1. $γ>0$⁠: risk-averse.

2. $γ=0$⁠: risk-neutral.

3. $γ<0$⁠: risk-loving.

The function is strictly increasing and concave for $γ>0$⁠, ensuring higher wealth leads to higher utility but with diminishing marginal returns. CRRA preferences are consistent with equilibrium results in the B-S model (Rubinstein, 1976; Brennan, 1979). By contrast, constant absolute risk aversion (CARA) utility is tractable in additive models, it assumes wealth-invariant risk aversion, conflicting with empirical evidence. Epstein–Zin preferences offer greater flexibility by separating risk aversion from intertemporal substitution but introduce significant computational complexity and are beyond the scope of this discrete-time setting. We acknowledge that CRRA does not fully disentangle time and risk preferences; its balance of realism, tractability and empirical validation makes it well-suited for our model. Indeed, CRRA preferences have been successfully applied to dynamic asset–liability management problems under model uncertainty, confirming its tractability and robustness in intertemporal decision-making contexts (Zhao, 2021). Future research may explore Epstein–Zin utility to generalize this framework in dynamic contexts.

Table 1 summarizes empirical estimates of the values of γ reported in the literature.

Table 1 presents a range of empirical RRA coefficients, some outside the domain adopted in our model. Several studies infer investors' risk preferences by contrasting observed market outcomes with risk-neutral benchmarks. For example, Bliss and Panigirtzoglou (2004) derived an adjusted risk-neutral probability density function incorporating a utility function to estimate the RRA from the FTSE 100 and S&P 500 option indices. Similarly, Ait-Sahalia and Lo (2000) introduced the Economic Value-at-Risk (E-VaR), a non-parametric risk measure based on state-price densities that better captures the market risk profile compared to traditional VaR.

While these contributions offer valuable insights, our model constrains the RRA coefficient (γ) to the interval −1 ≤ γ < 1, ensuring that the CRRA utility remains strictly concave, well-defined and computationally stable. We explicitly exclude the boundary case γ = 1, as it corresponds to a logarithmic utility function, a distinct model with different analytical properties that is not compatible with our current framework. Values outside this interval may result in excessive curvature or undefined CE, compromising both numerical consistency and interpretability. While Table 1 reports broader empirical estimates for completeness, all simulations and calibrations in this study are restricted within the specified domain. This modeling choice aligns with critiques of the risk-neutral valuation paradigm, such as those advanced by Friend (1977), reinforcing the behavioral and structural coherence of our approach.

Pratt (1964) established the foundations for analyzing risk aversion by introducing concepts such as risk premium, decreasing risk aversion and their proportional relation to total assets. Building on this, Fama (1970) showed that under uncertainty, even risk-averse investors maximize lifetime utility, and that nonstationary investment–consumption processes can be represented through non-stochastic single-period utility functions. Extending the utility framework, Richard (1975) proposed a multidimensional definition of RRA, highlighting consumption as a key determinant.

In parallel, Black and Scholes (1973) and Merton (1973) introduced asset pricing methods under a martingale framework, grounded in no-arbitrage and the uniqueness of the risk-neutral measure in incomplete markets. Thijssen (2008) illustrated this within a two-period economy, showing that future prices can be expressed as equivalent martingales and valued relative to present consumption.

Empirical evidence was provided by Brenner (2015), who analyzed 65,000 RRA coefficients from 7,000 U.S. executives (1996–2008) and found consistently moderate risk aversion (mean ≈ 1; median = 3). His framework-linked risk attitudes to firm characteristics, stock-option exercise behavior, liquidity needs and psychological effects such as overconfidence.

Ellersgaard and Tegnér (2018) integrated these theoretical and empirical strands by applying a martingale approach under the Heston (1993) model. They derived optimal CRRA-based consumption–wealth strategies in a stochastic volatility setting, empirically validated with both stocks and derivatives, including options.

In this section, we propose a primary framework based on the traditional differential stochastic equation, deriving new trend and diffusion terms that incorporate the CRRA utility function. Assuming that the underlying stochastic process follows a GBM, we employ the CRRA utility function to model the risk preferences of an economic agent, capturing both risk-aversion and risk-seeking behaviors.

In a risk-neutral setting, two critical conditions are essential for option valuation. First, estimated future cash flows are determined by discounting their expected values at the risk-free interest rate. A single neutral probability measure $(P*)$ must be derived from the subjective probability $(P)$ to prevent arbitrage on the discounted expected value. This ensures that the discounted price at an interest rate is a martingale and that the discounted expected value under this rate, according to the new probability $(P*)$⁠, does not present arbitrage opportunities. The probability measure used in a risk-neutral world is the martingale equivalent measure (MEM).

Consider the following traditional linear homogeneous differential equation ^[1]:

That describes the behavior of an specific underlying asset in a specific market, where $S0=S$⁠; $μ$ is a constant that denotes the asset's average rate return, $σ>0$ is the annual volatility, and ${Bt}t≥0$ is a standard one-dimensional Brownian motion (SODBM) define in the probability $(Ω,F,{Ft}t≥0,P)$⁠.

Consider the discounted pay based on the risk-free rate $r$ as $St*=e−rtSt$⁠. Then, apply Ito's lemma.

Allow us to get:

This process is supported by the theorem introduced by Heston (1993) and the application of the Cameron–Martin–Girsanov theorem as well as Girsanov's theorem (Mao, 1997). Then,

It is clear that $ϕ(t)∈L2([0,T];R)$⁠. Then, we obtain

Substituting (3) into (4)

Substituting (5) into (2), we got

It is worth to mention that $St*$ is Martingale respect $Ft$ to $0≤τ≤t$⁠. Now, replace (5) in (1).

Note that the previous equation corresponds to the behavior of an underlying price in a risk-neutral world.

Consider the traditional neutralized linear homogeneous differential equation (6) which describes the behavior of an underlying asset on a neutral-risk valuation, where $S0=S$⁠; $r$ is the risk-free rate, $σ>0$ is the annual volatility and ${Bt*}t≥0$ is a SODBM defined in the probability space $(Ω,F,{Ft}t≥0,P*)$⁠.

Assume the CRRA utility function $u(x)=x(1−γ)(1−γ)$⁠; with $γ≠1$ and $−1≤γ<1$⁠, considering this transformation $Yt=St(1−γ)(1−γ)$ and applying Ito's formula to $Yt$⁠, we obtain ^[2]:

with $α=(1−γ)(r−γ2σ2)$⁠; $θ=σ(1−γ)$⁠.

We recognize that incorporating risk aversion into a risk-neutral valuation framework may initially seem contradictory. However, our model resolves this by clearly distinguishing two separate valuation stages.

First, future stochastic payoffs are evaluated in a “risky world” or “risk-adjusted world” (Y), using the CRRA utility function to explicitly incorporate investor risk preferences. This process produces a risk-adjusted rate $α$⁠, which directly reflects the investor's degree of risk aversion. When the risk aversion coefficient $γ=0$⁠, $α$ coincides exactly with the risk-free rate, thereby recovering the standard risk-neutral benchmark.

Second, the adjusted payoffs are discounted to present value in a “risk-neutral world” (S), using the risk-free rate. Therefore, our approach accounts for risk in the payoffs themselves rather than in the discounting process, as summarized by the following formula:

here $r*$ represents the discount rate and $Φ(S)$ the final payoff.

This two-step valuation approach resolves the apparent contradiction by clearly separating the incorporation of risk preferences from the no-arbitration discounting condition. As a result, our model can be viewed as a generalization of traditional risk-neutral methods (e.g. B-S, CRR and LSM methods), which discount the expected weighted average payoffs directly under the risk-neutral measure (Björk, 2009).

To build a multiplicative tree that allows us to describe the discrete behavior of the underlying asset, we must determine the first two moments of the ordinary linear stochastic differential equation to obtain the transition probability, which we present next. These expressions have been widely employed to calibrate recombining lattice models (Lari-Lavassani et al., 2001; Marín-Sánchez, 2010), which we adopt in this work ^[3].

The match between the first and second moments for both processes allow us to obtain a specific recombination given by

where $AYti,Δt:=EYt∣Yti⌋ti+ΔtYji;BYti,Δt:=EYt2∣Yti⌋ti+ΔtYji2;CYti,Δt:=1+BYti,Δt2AYti,Δt$

Proof. Using the definition of $Ed[·|·]$⁠, the matching between the discrete and continuous time processes is

Because $Yj+1(i+1)=Yj(i)uj(i)$ and $Yj(i+1)=Yj(i)dj(i)$ the previous equations can be reduced to

The previous equations can be rewritten as

From Equation (12) we obtain:

Substituting Equation (14) into Equation (13) we obtain:

Without loss of generality and to eliminate a degree of freedom, we assume $uj(i)dj(i)=1$ therefore, from Equation (15) we obtain

If we assume that, based on $C(Yti,Δt)$ we can derive an expression to determine $uj(i)$⁠, we obtain the following expression:

And finally

From Proposition 2 and Equation (9), we find

Similarly, using Equation (11), we find

Hence

So,

Substituting Equations (18) and (19) on (16), and taking the positive solution of the square root, we obtain

Finally, we calculate $pj(i)$ from Equations (17) and (20), which are substituted for Equation (14) as follows

Consider a stochastic process $St$ modeled discreetly over time. In each time period $i$⁠, the possible states of the process are given by a finite number of vectors $S(i)=(Sj(i))$⁠, where $j$ is a finite indexed set $J(i)$⁠. Just as the time stream runs from $i$ to period $i+1$⁠, discrete node $Sj(i)$ can become part of any $Sj(i+1),j∈J(i+1)$⁠.

The discrete process determines a transition probability matrix $P(i)=Pj,j(i)$⁠, and thus the element $Pj,j(i)$ represents the probability that $Sj(i)$ moves to $Sj(i+1)$⁠. The number of rows in $P(i)$ is the cardinal of $J(i)$ and the number of columns is the cardinal of $J(i+1)$⁠, such that all $j∈J(i)$ and $∑j∈J(i+1)Pj,j(i)=1$⁠.

The transition probabilities are given for each i by a vector $Pj(i),j∈J(i)$ such that $Pi,j(i)$ is defined by $Pj(i)$ if $J=j+1,1+Pj(i)$ if $J=j$⁠, and zero otherwise.

In this multiplicative process, the size of the upward and downward jumps is defined by $uj(i)$ and $dj(i)$⁠, respectively, so that $Sj+1(i+1)=Sj(i)uj(i)$ and $Sj(i+1)=Sj(i)dj(i)$ such that $uj(i)dj+1(i+1)=dj(i)uj(i+1)$⁠; obtaining the recombination given by

Note that in the special case where $uj(i)=u$ and $dj(i)=d$ are constants, this recombination takes the form $Sj(i)=S0uj−1di−j+1$⁠.

In this section, we provide a brief overview of the MBRT with CRRA (MBRT–CRRA) method for option valuation, incorporating a CRRA-type utility function that corresponds to the subjective probability $P$⁠. This probability is defined in the real-world measure, accounting for arbitrage and the absence of risk neutrality. Finally, we present the necessary transformation that enables evaluation under the risk-neutral measure.

Assume that the financial option of an asset with an initial value of $Y0=S0(1−γ)1−γ$ and exercise price of $K$⁠, and $T$ is divided into $N$ subintervals, each with a length of $Δt$⁠. We define $fj(i)$ as the value of the option in the $(i,j)$ node. Based on Marín-Sánchez (2010), the price of the asset, expressed as a general function, has binomial recombination in the $(i,j)$ node, which can be represented by the following equation:

with $Y0=Y1(1);i=0,1,…,N;j=1,2,…,i.$ It is important to note that, in this case, both $u$ and $d$ are constants Lari-Lavassani et al. (2001), so Equation (22) is summarized as follows:

This is followed by an algebraic approach derived from the proposed methodology to assess the basic financial options.

In case of a European call option, the value on the maturity date is indicated by $(max(St−K,0))1−γ1−γ$⁠, where $St=((1−γ)Yt)11−γ$⁠, recall that $γ≠1$ and $−1≤γ<1$⁠, so

whereas the value at each node is given by

The same analysis can be applied to valuing European put options. In the case of an American call option, the value at maturity is estimated as the value of a European call option, while the discounted value is defined by

Similarly, we can represent the value of the maturity date for an American put option. In the case of a Bermudan call option, the value at maturity is estimated in the same way as for an American or European option, but in some switching periods, $fj(i)$ changes from Equation (24) to Equation (23).

To ensure that the option valuation does not yield incorrect or spurious solutions at the nodes of the binomial tree, it is essential to verify the numerical properties of the proposed schemes. For instance, we verify whether the probability of transition lies within the interval $[0,1]$⁠, whether the binomial recombination representing prices remain non-negative, and whether the valuation scheme is positive a monotonically increasing in the time direction. Once these conditions are identified, we establish specific constraints for the scheme that adhere to ideal properties, which are maintained through a series of propositions proven in the next section.

Suppose a Bermudan call option at different maturities, where the transition between the European and American payoff occurs at $T/2$⁠. Let us assume that the dynamic behavior of the of a given financial asset follow a GBM, as described in Equation (8), with the initial conditions $S0=10$⁠, $K=10$⁠, $γ=0$⁠, $r=0.05$⁠, and annualized volatility of $σ=0.30$⁠. These values allow us to calculate the parameters $α$ and $θ$⁠. In this framework, we consider that investors exhibit either risk-averse or risk-loving behaviors, which are captured by the RRA coefficients $γ=0.3$ and $γ=−0.3$ respectively. For reference, we also present the classical risk-neutral case. To evaluate the performance of our approach, we compare the experimental numerical results obtained using the MBRT–CRRA and LSM Monte Carlo Simulation with CRRA (LSM–CRRA ^[4]) methods. In the LSM–CRRA method, the option value is estimated based on 10,000 simulated trajectories and 1,000 repetitions.

Table 2 shows that option values rise with longer maturities, and both methods yield consistent results, supporting the validity of our framework. Importantly, risk preferences alter premiums: under risk aversion, values are lower than in the risk-neutral case, reflecting the tendency of investors to exercise early to avoid uncertainty. By contrast, risk-loving investors hold positions longer, increasing counterparty exposure and leading to higher premiums relative to the neutral benchmark.

Using the initial conditions of $S0=10$⁠, $N=100$⁠, $σ=0.3$ and $r=0.05$ in Equation (8), Figure 1a illustrates that the value of $α$ tends toward zero as coefficient $γ$ approaches +1, but remains positive if $γ$ differs from +1. Under the B-S valuation, this is exactly 5%, consistent with the assumed risk-free interest rate and corresponds to $α$ is zero meaning that the investor exhibits risk neutrality.

Figure 1b further illustrates that changes in the RRA coefficient directly affect profit price fluctuations. For positive values up to 0.57, profits decline and may even turn negative; beyond this threshold, they rise again, implying that the model should be constrained within this range to maintain economic coherence.

By contrast, negative values of the coefficient (risk-seeking behavior) lead to steadily increasing profits as risk preference intensifies. It is important to note that when $γ=1$⁠, a discontinuity arises. This issue can be easily resolved using L’Hôpital's rule by transforming the function into a logarithmic form.

Figures 2a and b illustrate the evolution of both $St$ and $Yt$ over time, assuming variations in the parameter $α$⁠. Figure 2a demonstrates that the price of the underlying asset and its corresponding utility exhibit a similar patterns. The key distinction is that when the utility coefficient is negative, the value of the utility is higher; conversely, when the RRA coefficient is positive, the utility decreases. Nonetheless, the relative changes remain proportional.

Figure 3a presents the results obtained using the LSM–CRRA and MBRT–CRRA methods, illustrating the convergence between our proposed approach and the B-S methodology for valuing a European call option. The valuation is performed under the initial conditions of $S0=10$⁠, $K=10$⁠, $N=100$⁠, $σ=0.3$⁠, $r=0.05$⁠, $T=1$ years, $γ=0$ and 50 repetitions of 10,000 trajectories.

The MBRT–CRRA method converges rapidly to the exact solution within 40–50 periods, while the LSM–CRRA approximation oscillates below the benchmark and stabilizes after 60 periods.

The numerical results show that the MBRT–CRRA provides an adequate price compared to the LSM–CRRA method, as the difference between the two methods is not significant. Additionally, in terms of time efficiency, the MBRT–CRRA shows an outstanding performance with reductions between 96 and 99% in price estimation. In terms of computational performance, the MBRT–CRRA method requires approximately one second to compute option price with N = 300, by contrast, the LSM–CRRA method takes about three seconds for a single run (10,000 trajectories) and, given its reliance on simulation, must be repeated multiple times – at least 50 repetitions – to obtain stable results.

The figure also shows that LSM–CRRA provides a high-accuracy estimate, with confidence bands indicating that deviations from the exact value remain within 5%.

Figure 3b, shows that MBRT–CRRA and LSM–CRRA converge across different values of γ. For positive coefficients (e.g. +0.8) convergence is reached quickly due to the stability of LSM–CRRA relative to the dispersion of the MBRT–CRRA. When $γ$ is negative, the convergence between the methods occurs near the end of the proposed interval. In this scenario, the MBRT–CRRA method proves to be more stable while the LSM–CRRA tends to exhibit higher volatility, particularly with a coefficient of 0.8. In both cases, the exact value cannot be derived using B-S, since it lacks risk-adjustment terms, underscoring the motivation for our framework and the scope for deriving an exact differential equation in future research.

The primary objective of this study is to present a numerical framework that incorporates risk aversion into option valuation. To this end, we conducted four sensitivity analyses to examine how the MBRT–CRRA and LSM–CRRA methods capture investor risk preferences and their impact on option premiums. These analyses were performed with $S0=10$⁠, $K=10$⁠, $r=0.05$⁠, annualized volatility of $σ=0.30$⁠, $T=1$⁠, $N=300$⁠, and using the LSM–CRRA method with four regressors.

Figure 4 clearly illustrates how the variation in annualized volatility $(σ)$ and the $γ$ coefficient affect option prices across European, American and Bermudan (call and put options) under the LSM–CRRA and MBRT–CRRA frameworks. The sensitivity analysis demonstrates a pronounced asymmetry: when investors are risk-loving, they consistently assign higher values to options as volatility increases. Conversely, for risk-averse investors, option prices initially decrease with low volatility and only begin to rise again at higher volatility levels. For example, at $γ=+0.1$⁠, option values initially decline within the 1%–3% volatility range but subsequently increase as volatility rises beyond this threshold. As $γ$ approaches higher levels of risk aversion (e.g. $γ=+0.7$⁠), this recovery of option values is delayed, indicating a stronger aversion to volatility.