Numerical solution of a nonlinear PDE model for pricing Renewable Energy Certificates (RECs)

In this article we present a valuation method for Renewable Energy Certiﬁcates (RECs) or green certiﬁcates. For this purpose, we propose a non-linear PDE model with two stochastic factors: the accumulated green certiﬁcates sold by an authorized generator and the natural logarithm of the renewable electricity generation rate. One novelty of the work comes from the numerical treatment of the non-linear convective term in the PDE. Thus, we use the Bermdez-Moreno algorithm to deal with this non-linear term. This duality algorithm is based on the Yosida regularization of non-linear maximal monotone operators. In order to solve the obtained linearized problem, we use numerical methods based on semi-Lagrangian schemes in the direction without diffusion while we consider an implicit second order ﬁnite differences scheme in the direction with diffusion term. Finally, we show illustrative results of the performance of the proposed model and numerical meth- ods that have been implemented.


Introduction
In recent years, several governments have encouraged environmental policies for promoting renewable energy sources and have established targets for renewable energy growth. Nevertheless, renewable energy technologies often require a large investment, hence alternative market tools are required to achieve these targets. One tool to implement these policies and to incentivize investing in renewables is the use of renewable energy certificates or RECs (often called green certificates or GCs in Europe). When renewable energy generators meet certain criteria, they receive one certificate for a specific unit, typically 1 MWh of renewable electricity produced. This REC can be sold to a Load Serving Entity (LSE) that is subject to the annual requirement of some percentage of electricity procured from renewables. If this quota obligation is not met by the LSE, a non-compliance penalty applies. This penalty is called the Alternative Compliance Payment (ACP).
Many countries or states have adopted renewable portfolio standards (RPSs) and trading of RECs. These instruments can also be linked to the generation of a particular type of energy, as is the case of the New Jersey (NJ) market for SRECs (solar renewable energy certificates) presented in [13] . Similar markets exist around the world (eg, Italy, UK, Sweden, Norway, US or Australia) and are considered an important alternative to implementing other environmental policies such as taxes or regulatory limits, thus helping to address climate change, especially in countries with an absence of carbon emissions markets, such as the US.
The academic treatment of REC markets is more recent than the study of cap-and-trade schemes for carbon emissions (see [9 , 12 , 16 , 17 , 22 , 26 , 27] ). Although related, there are several differences between carbon and REC markets: supply and demand in opposite roles, different inter-temporal connecting mechanisms (such as withdrawal or unlimited banking) or external and underlying factors (such as fuel prices, power demand or renewable energy generation). Nonetheless, related approaches exist for the mathematical modelling of the valuation problem of green certificates. Some models of REC markets are obtained by using a stochastic dynamic setting, thus replicating the price volatility (see [21] and [18] ).
In the present work, we address the valuation of green certificates as the solution to a coupled system of forwardbackward stochastic differential equations (FBSDE). This type of equations have been used for the pricing of financial instruments in emission markets in [10] and [17] . More precisely, we consider a forward stochastic differential equation (SDE) for the renewable generation rate and another one for the accumulated green certificates, and a backward stochastic differential equation (BSDE) for the green certificate price. This backward component of the FBSDE can be formulated in terms of a semilinear partial differential equation (PDE). As the analytical expression of a solution for the PDE is not available, in the present article we propose a set of numerical techniques to solve this non-linear equation. More precisely, the difficulty associated to the non-linear convection term has been solved by means of a duality method, which is known as Bermdez Moreno algorithm (see [6] ). This algorithm is based on the Yosida approximation of non-linear maximal monotone operators. In [2] this method has been applied to solve a non-linear diffusion problem, while in [7] it has been used for the multivalued maximal monotone operator arising in the solution of complementarity problems related to Asian options with early exercise opportunity, also referred as Amerasian options. At each step of this algorithm, we need to solve a degenerated linear PDE problem. For this purpose, we use a characteristics scheme to discretize the material derivative operator in one of the spatial directions combined with the use of a second order implicit finite differences scheme in the other spatial direction.
This article is organized as follows. In Section 2 , we introduce the mathematical modeling corresponding to the two stochastic factors. In Section 3 , we describe the numerical methods, including the treatment of the non-linear convective term, the analysis of boundary conditions and the full discretization of the PDE. Section 4 shows some illustrative numerical results. Finally, some conclusions are presented in Section 5 .

Mathematical modelling
In order to model the price of a REC we assume that it is given by a stochastic process, which is denoted at time t by P t and depends on two stochastic factors. More precisely, these factors are the renewable generation rate and the number of accumulated green certificates, which are denoted by G t and B t , respectively. We assume that the values of these two factors are known at the initial time t = t 0 .
Specifically, for t ∈ [ t 0 , T ] , the generation rate G t is most influenced by weather patterns and by the construction of new renewable capacity. It is therefore natural to assume that the dynamics of G t is given by where ˜ G t is a Ornstein-Uhlenbeck (OU) process which satisfies the following stochastic differential equation (SDE): where μ g (t, ˜ G t , P t ) and σ g are the drift and the volatility, respectively, and W 0 t is an F t -adapted Q -Brownian motion, where Q denotes the probability measure. In particular, we assume that the drift is linear with respect to P and given by the expression where α g is the mean reversion speed of the process and β g is the parameter which controls the level of feedback from the price of the certificate. This feedback parameter captures the elasticity of supply to price, meaning the tendency of new renewable generation to be installed at times when certificate prices are high. Moreover, the seasonality is represented by the deterministic function f, which depends on time t. This function f is a combination of cosine and sine functions to represent the influence of weather on prices. For this purpose, there are more options to represent these seasonal patterns, such as the time series used in the particular case of Swedish-Norwegian market (see [4] ).
Similarly to some methodologies used in the pricing of Asian options, the number of accumulated green certificates at time t, B t , satisfies the following equation: Moreover, as accumulation is measured from the beginning of the compliance period, t 0 , we assume that B t 0 = 0 . Note that since the process B t represents an accumulated quantity, it turns out to be positive and non-decreasing.
The main objective of the model is to characterize the price of the renewable energy certificate, P t , at time t. In our formulation in terms of a PDE problem, we assume the existence of a function P, such that P t = P (t, B t , G t ) . Once the function P is obtained as the solution of the PDE problem, we can compute the value of P t for given values of t, G t and B t .

Single period case
We denote by T the maturity of the certificate and by γ the number of life years of the certificate (i.e. the number years since issuance for which the REC is valid to submit for compliance). Then, we initialize the life of the certificate at t 0 = T − γ . In the single period case (i.e. assuming one year with no intermediary compliance before T ), we take γ = 1 . The value of the REC at time t = t 0 is unknown. However, its value at maturity, P T , is given by the terminal condition: where ψ : R + −→ R + is a bounded measurable and decreasing function. More precisely, it is given by ψ ( where π T is the penalty amount at time T and R T is the requirement at time T . Since the discounted price of the REC is a martingale under Q (a 'no arbitrage' condition), the price of the certificate at time t is equal to the discounted value of the conditional expectation of its terminal value, i.e.
where r is the constant risk free interest rate. The previous expression implies that the process P t is bounded, taking values in [0 , π T ] . Moreover, since the filtration F t is generated by the Brownian motion W 0 t , the price of the renewable energy certificate can be represented as an Ito integral with respect to W 0 t as follows for some F t -adapted square integrable process Z 0 t . Next, by using Eqs. (2) , (4) and (7) , the price process of a renewable energy certificate P t in a single period can be described by the solution of the following FBSDE in the time interval [ t 0 , T ] : Note that there are two kinds of stochastic differential equations in (8) depending on the direction of time: the first two equations are forward ones while the last is a backward equation. Analogous equations have been considered for carbon emission prices in [17] .
Assuming the existence of a solution for (8) and where P t is a traded asset with a drift equal to risk neutral rate under the risk neutral measure (as indicated in the third equation of (8) ), we can use Ito's formula for a process depending on the two Ito processes B t and ˜ G t (see [19] , for example) to obtain: Therefore, identifying the drift coefficient of the previous expression and the corresponding one of the third equation in (8) , the function P = P (t, B, ˜ G ) satisfies the following non-linear PDE: Taking into account that the value of the REC at maturity is given by (5) , then PDE (9) jointly with the final condition defines the final value problem for the single period case.

Multiple periods case
The previously presented arguments for the single period case can be extended to an arbitrary number of compliance periods. In the multiple period case, the price of the certificate at maturity T is equal to the payoff (5) , while a jump condition must be applied at each compliance date T i , i = 1 , . . . , γ − 1 . Thus, the value of the certificate at the compliance date T i is given by where ψ (·) = π i 1 [0 ,R i ) (·) and R i is the requirement at time T i . In this case, a sequence of linked final value problems is defined by Eqs. (9) and (11) . The existence and uniqueness of solution for the final value non-linear PDE problem (9) - (10) is an open question (analogously for the final value problems (9) -(11) ), that could be addressed in a future research. Although this issue is not addressed this article, first we note that the FBSDE (8) is non-standard and similar to the one appearing in emission markets in [17] , the existence and uniqueness of solution of which is rigorously analyzed in [25] . So, the ideas in [25] could be used to address existence and uniqueness of solution of (8) . Concerning the existence and uniqueness for the PDE problem (9) -(10) , the governing PDE is parabolic degenerate in one direction and semilinear. As indicated in [25] for the case of emission certificates, the connection between the FBSDE and the PDE is based on an extension of the Feynman-Kac formula that appears in [23] and that requires the regularity of the classical solution of the PDE, which we are not able to prove yet.

Numerical methods
Once the final value problems associated to the single period and multiple period cases have been posed in the previous section, since there are no analytical expressions for their respective solutions, we propose a set of numerical methods to approximate them.
For this purpose, first note that the PDE (9) is initially posed in a spatial unbounded domain, so that a truncation of the domain to a bounded one is required for the numerical solution. This domain truncation is very common in option pricing problems, for example, and requires an appropriate selection of the truncation boundaries and the corresponding boundary conditions to obtain a proper approximation of the solution in the region of actual financial interest [20] . Secondly, as the drift term (3) depends linearly on the unknown P, the corresponding first order (convective) term introduces a nonlinear aspect that needs to be solved. In the present work we will take advantage of the fact that the non-linear term can be written in terms of a maximal monotone operator to apply a duality method. A third aspect is related to the lack of second order derivative with respect to variable B, which makes the second order PDE degenerate. In order to overcome the difficulties associated to this last aspect, we propose a semi-Lagrangian numerical scheme to jointly discretize the terms of time derivative and the first order derivative with respect to B .
Clearly, in view of these important specific characteristics of the PDE (9) suitable and efficient numerical techniques need to be applied.

Treatment of the non-linear convective term
As previously pointed out, the PDE problem (9) includes a non-linear convective term. One possibility to deal with this non-linearity is based on the Bermdez-Moreno algorithm, which involves the Yosida regularization of non-linear maximal monotone operators (see [6] ). In [2] , the here proposed techniques have been applied to a one-dimensional problem with non-linear diffusion term instead of the non-linear convection one. Next, in [1] the convergence of the method has been analysed and in [3] extended to two spatial dimensions.
Thus, following the idea used in [2] , let us first introduce the maximal monotone operator m, defined by so that Therefore, Eq. (9) can be written in the equivalent form:

∂P ∂t
Following the duality technique developed by Bermúdez and Moreno in [6] , for the constant parameter ω > 0 , we introduce the new additional unknown where I denotes the identity operator. So, we have 1 2 Next, by using the Bermúdez-Moreno lemma from [6] , also applied in [3] , for λω < 1 , we can obtain the equivalence where m ω λ denotes the Yosida approximation of the operator m ω = m − ωI with parameter λ. This Yosida approximation is is the resolvent operator of m ω , which is defined for λω < 1 and it is a monotone Lipschitz function with constant (1 − λω) −1 . It is easy to prove the equivalence (17) . For this purpose, starting from its right hand side of (17) , first we have: which is equivalent to or equivalently: this last identity being the left hand side of (17) . For details about Yosida approximation of maximal monotone operators, we address the reader to [8] . In [1] , the convergence of the method for a model with the same nonlinear term in the diffusion part of an elliptic operator in one dimension has been analyzed, obtaining that the optimal choice for convergence comes from condition 2 λω = 1 . Although we have not addressed this theoretical analysis for the here treated problem, we consider also this relation between both parameters. Also this choice is theoretically obtained in [6] for some elliptic variational inequalities formulated in terms of a multivalued subdifferential operator. Under this choice, the expression of Yosida approximation can be computed and is given by [2] : Next, if we introduce the linear differential operator and take into account the expressions (15) and (16) , then Eq. (14) can be rewritten in the equivalent form: Moreover, from the equivalence stated in (17) , Eq. (19) is coupled with the following non-linear equation: In order to solve the non-linear system given by (19) - (20) , we propose a fixed point iteration which mainly starts with an initial value of θ to solve the linear PDE (19) . Then, we replace the obtained P and the last computed value of θ in the right hand side of (20) to update the value of θ and start solving again the PDE (19) in the next step. This fixed point algorithm will be explicitly written in a forthcoming section once the discretized problem has been introduced.

Localization and analysis of boundary conditions
In order to apply the numerical discretization using finite differences to the PDE problem, it is necessary to consider a bounded computational domain. Thus, for a given value of θ we approximate the linear PDE problem (19) In order to establish the boundaries of the truncated domain which require boundary conditions to be imposed, we follow the results in [24] and introduce the following notation and we use the notation * = Next, we denote by n = (n 0 , n 1 , n 2 ) the normal vector to * pointing inwards * . Let us define the following subsets of * : a ij n i n j = 0 , ∂a ij ∂y j n i < 0 .
Following [24] , we need to impose boundary conditions at 1 ∪ 2 . Thus, for (21) , we conclude: Hence, we impose the following homogeneous Neumann boundary conditions at the spatial boundaries: jointly with the final condition at the boundary y 0 = T in the single period case. In addition, in the multiple period case the condition (11) is applied at the boundaries y 0 = T i , which represent each compliance date.

Discretization of the PDE
In order to choose an appropriate time discretization scheme for the PDE (19) , we note that the linear differential operator (18) is degenerate. Therefore, PDE (19) can be considered as a limit case of a convection dominated PDE, especially in the direction without diffusion term. Therefore, in order to avoid spurious oscillations related to the use of standard finite differences schemes, we propose a suitable version of the method of characteristics (see [5] or [14] , for example). Note that PDE (19) turns out to be similar to the one arising in the pricing of Asian options with continuous arithmetic averaging. Thus, we follow the idea first proposed in [15] for this kind of problem, which consists of choosing a semi-Lagrangian method (also referred as the characteristics method) in the direction without diffusion combined with a Crank-Nicolson finite differences scheme in the direction with diffusion. This method has been also used for solving PDE models for the valuation of business companies in [11] . Note that in our coupled non-linear problem (19) -(20) we apply the proposed time discretization scheme to the step of solving Eq. (19) .
For the time discretization, we first consider the change of time variable τ = T − t, where τ represents the time to maturity. Therefore, Eq. (19) can be equivalently written in the domain ˜ = (0 , γ ) × (0 , 1) × (0 , 1) as follows where the involved differential operators are given by Note that D P D τ represents the material derivative of P in the direction ˆ B associated to the one-dimensional velocity field , which does not depend on B . Moreover, A denotes the second order convection-diffusion-reaction differential operator in the direction ˆ G . Note that this splitting of the differential operator governing equation (19) is the departure point of the proposed time discretization scheme.
First, we use a characteristics scheme to discretize in time the term associated to the material derivative. This semi-Lagrangian scheme is based on a finite difference discretization of the time derivative along the characteristic lines (see [5 , 14 , 28] ). For this purpose, we introduce N T > 0 and a time step τ = γ /N T for considering a uniform time mesh with nodes given by τ n = n τ for n = 0 , 1 , . . . , N T . At each time step we consider the initial value ODE problem satisfied by the trajectory associated to the velocity field v through the point (τ n +1 , ˆ Note that the solution of this ODE problem is given by In order to build the finite differences approximation of the material derivative along the characteristics, we introduce χ n = χ (τ n ) , which is given by and represents the position at time τ n of the point placed at ˆ B , ˆ G at time τ n +1 and moving according to the velocity field v .
Next, we introduce the approximation for the material derivative: Secondly, by using a Crank-Nicolson scheme ( ˆ θ = 0 . 5 in the so called ˆ θ -method) for the second order differential term A P in Eq. (25) , we obtain: At each time step, the evaluation of the term P n • χ n in (27) at the quadrature nodes is approximated by using a biquadratic interpolation formula from the values of P n at the mesh nodes.
Note that at each time step, Eq. (27) is coupled with the following non-linear relation between P n +1 and θ n +1 : Next, we propose a fixed point algorithm to approximate the solution of the non-linear problem (27) - (28) . This fixed point algorithm mainly consists of solving Eq. (27) to obtain P n +1 for a previously computed value of θ n +1 , and next updating θ n +1 according to (28) with the more recent values of P n +1 and θ n +1 . Thus, the algorithm can be sketched as follows: 1. Let P 0 and θ 0 be initialized (for example θ 0 = 1 ).
• For a given θ n +1 ,k , we obtain P n +1 ,k +1 by solving jointly with the boundary conditions.
• We check the stopping test (c) If the stopping test is satisfied then go to 2, otherwise go to (b).
In order to describe the solution of the fully discretized problem, let us introduce the notation (τ n , ˆ to represent a generic node of the uniform finite differences time-space mesh with time step τ and spatial steps ˆ B and ˆ G , for indexes n = 0 , 1 , . . . , N T , i = 0 , 1 , . . . , N ˆ  B and j = 0 , 1 , . . . , N ˆ G . At each fixed point iteration, the full discretization of problem (29) can be written as follows: where ˆ θ = 0 . 5 for the Crank-Nicolson time discretization, and P l,m r,s ≈ P m (τ l , ˆ By taking into account the previous expression of the fully discretized problem, we have to solve a linear system with (N ˆ B − 1) × (N ˆ G − 1) unknowns at each time step. Moreover, if we order the finite differences mesh nodes in lexicographical order, the resulting matrix is block diagonal with N ˆ B − 1 blocks of tridiagonal matrices of order N ˆ G − 1 each. So, by applying the classical Thomas algorithm for tridiagonal matrices, each N ˆ B − 1 linear system can be efficiently solved. Thus, at each time step and for each value of i = 1 , . . . , N ˆ B − 1 , we have the following linear system: i,N ˆ G −1 is the appr oximation of the solution at the finite differ ences mesh nodes with coordinate ˆ B = ˆ B i , and the matrix C( ˆ G ) is given by Moreover, for j = 1 , . . . , N ˆ G − 1 , the jth component of the second member vector b n i of the linear system is given by

Numerical examples
In this section we present some numerical results to illustrate the performance of the proposed numerical method, as well as to discuss some quantitative and qualitative results for a real problem.

Academic test
As a sanity check of the code and numerical methods, in the first example we show an academic test with known analytical solution. For this purpose, we consider the following non homogeneous non-linear PDE: where the differential operator L 1 is defined by (9) and h is given by   In this first academic test we do not include the seasonality effect, so that we take f = 0 . Parameters in the PDE are collected in Table 1 and mostly taken from [13] .
For the duality method we consider the parameter ω = 2 and the stopping test = 10 −5 . In order to assess the performance of the proposed numerical methods, we consider a constant relationship between the time and spatial steps, i.e.
G . Next, we compute the discrete L ∞ relative error between the exact solution and the numerical approximation at the time step n as follows: Next, we consider the maximum of errors defined in (31) , Er r ∞ ( τ ) = max n (er r ∞ n ( τ )) . Moreover, the radius of the convergence is given by for the stepsize τ, and the empirical order of convergence is given by log 2 (R ) . Table 2 shows the errors, convergence ratio and empirical order of convergence with different time and spatial discretizations computed as in [15] .
Thus, taking into account the relative errors and the convergence ratio in Table 2 , we can conclude that a first order convergence is achieved.

Real case
In this example we analyze the evolution of the price of a real green certificate. For this purpose, we have used the real New Jersey market data presented for SREC markets, i.e., markets for solar renewable energy certificates, in [13] .
In this market, the energy year refers to the 12-month period ending on May 31. We assume that the maturity is T = 13 , i.e., May 31, 2013. For convenience, the energy year 2013 is defined as the time interval (12,13]. Thus, we consider the requirement schedule 2010-2013, i.e., γ = 3 , with initial year t = t 0 = T − γ = 10 and the first compliance date at t = 11 (at the end of the year 2011) which corresponds to τ = 2 . For the PDE parameters we consider those ones in Table 3 .
Moreover, the seasonality function f represents the influence of weather conditions and is chosen as follows: f (s ) = a 1 sin (4 π s ) + a 2 cos (4 π s ) + a 3 sin (2 π s ) + a 4 cos (2 π s ) , the parameters of which are provided in Table 4 . Table 5 shows the requirement, R i , and the penalty, i , values at the end of each year for i = 1 , 2 , 3 .
For the numerical methods, we start by choosing ˆ b = 7 × 10 5 and ḡ = ln (7 × 10 5 ) , so that ˆ g = 2 × ln (7 × 10 5 ) . By using these values we pose the PDE problem in the bounded domain ˆ = (0 , γ ) × (0 , 1) × (0 , 1) . Concerning the discretization parameters, we consider 100 time steps per month for the time discretization, i.e. τ = 1   Next, we show the computed results in the variables (t, B, G ) that are obtained with the model data and the parameters related to the numerical methods. Firstly, in Fig. 1 we show the price of the certificate eight months before maturity, that is at time t = T − 2 / 3 (or at τ = 2 / 3 , equivalently). We can observe that the price of the certificate takes values between zero and the penalty amount. When the accumulated supply (B) becomes very high, the price of course becomes very low, and when supply is very low, the price turns out to be very high. Moreover, when the number of accumulated certificates (B) and the renewable energy rate (G = exp ( ˜ G )) tend to zero, the price approaches to penalty value π . That is, for low values of both variables, B and G, the price of the certificate is almost equal to the penalty. Moreover, we also point out that the price of the certificate decreases when we increase the value of both state variables since, as it is indicated in [13] , for high values of the variables, supply of certificates is high so the market can achieve the requirement easily and an additional certificate will not be needed for compliance. Therefore, when only B tends to zero, the price tends towards the penalty unless high values of G offset this, as it is illustrated in Fig. 1 . When we represent the value of the certificate for a time closer to maturity, for instance at time t = T − 1 / 3 (i.e. four months before expiry date), as is shown in Fig. 2 , we observe that low values of generation rate are associated with prices equal to the penalty amount even for values of accumulated certificates nearer to requirement since time is running out before compliance. High values of generation rate are linked to prices equal to the penalty only for lower values of banked certificates. As expected, a high current production of renewable energy compensates for a low number of accumulated renewable energy certificates.
Finally, if we want to show the price of the certificate versus the accumulated renewable energy certificates for different times, we can create cross-sectional plots as the ones in Fig. 3 . In order to obtain these curves we have chosen appropriate large values of renewable energy rate. At maturity, t = T , the price of the certificate is equal to the penalty if the requirement is not met, otherwise the price is zero. Then, as we move backwards in time we can observe that the curves move to the left and take lower values, due to the diffusion of the final value. When we arrive at a compliance date there is a jump to the right and the price increases again due to the jump condition imposed. In such figure we show the value at the first compliance date as well. Note that at maturity and compliance dates there exists a discontinuity when the number of accumulated certificates meets the requirement, because of the indicator function multiplied by the penalty, as it is pointed out in [13] .
The CPU time is 223 seconds for this real case example. The CPU of the laptop we use is a Intel(R) Core(TM) i7-10750H at 2,6 GHz with 32 GB 2933 MHz DDR4 RAM. The implementation has been developed in Matlab.

Conclusions
Climate change policy and the growth of renewables have led to the creation of many new tradeable certificate markets for which price modelling and derivative pricing remains particularly challenging. Complicated market designs produce feedback mechanisms whereby higher certificate prices lead to higher incentives to build more renewables, ultimately increasing long-term certificate supply, and lowering expected prices. We have provided in this paper a new PDE formulation for the resulting equilibrium price, building on earlier work on REC or green certificate prices. As we have seen, an interesting PDE problem emerges for the certificate price in this setting, with a non-linear convective term in a PDE which is also degenerate. A set of numerical methods is proposed to handle these issues, and the performance of the approach is confirmed via a real world example, taken from the NJ SREC market. Such a model is particularly valuable for renewable energy investors to better understand their certificate price risk, as well as the potential benefits of a new investment. It may also be an informative tool to regulators and policy makers and they continue to improve market design to further incentivise global renewable energy growth.