We get the same LP form except for the set of inequality constraints:

$$ C x \leq D \Leftrightarrow \left( \begin{array}{ccc} I_n & -I_n & \mathbf{0}_{n,n_{\mathcal{S}ector}} \\ -I_n & -I_n & \mathbf{0}_{n,n_{\mathcal{S}ector}} \\ C_{\text{DTS}} & \mathbf{0}_{n_{\mathcal{S}ector},n} & -I_{n_{\mathcal{S}ector}} \\ -C_{\text{DTS}} & \mathbf{0}_{n_{\mathcal{S}ector},n} & -I_{n_{\mathcal{S}ector}} \\ \mathbf{CI}(t)^{\top} & \mathbf{0}_{1,n} & \mathbf{0}_{1,n_{\mathcal{S}ector}} \\ \mathbf{CM}(t)^{\top} & \mathbf{0}_{1,n} & \mathbf{0}_{1,n_{\mathcal{S}ector}} \\ -\mathbf{GI}(t)^{\top} & \mathbf{0}_{1,n} & \mathbf{0}_{1,n_{\mathcal{S}ector}} \end{array} \right) x \leq \left( \begin{array}{c} b \\ -b \\ \text{DTS}^{*} \\ -\text{DTS}^{*} \\ (1 - \mathbf{R}(t_0, t)) \mathbf{CI}(t_0, b(t_0)) \\ \mathbf{CM}^{*}(t) \\ -(1 + \mathbf{G}) \mathbf{GI}(t, b(t)) \end{array} \right) $$

and the upper bound: $x^+ = (1 \{ \mathcal{C}M(t) \leq \mathcal{C}M^+ \}, \infty \cdot \mathbf{1}_n, \infty \cdot \mathbf{1}_{n_{\mathcal{S}ector}})$.

**Example 62** We consider Example 60. The carbon momentum values are equal to -3.1%, -1.2%, -5.8%, -1.4%, +7.4%, -2.6%, +1.2%, and -8.0%. We measure the green intensity by the green revenue share. Its values are equal to 10.2%, 45.3%, 7.5%, 0%, 0%, 35.6%, 17.8% and 3.0%. The net-zero investment policy imposes to follow the CTB decarbonization pathway with a self-decarbonization of 2%, and to improve the green intensity of the benchmark by 100%.

The solutions of the LP optimization problem are shown in Table 11.25. Note that again there is no solution for some dates. In addition, compared to the net-zero equity portfolio exercise, we have changed the self-decarbonization target $\mathcal{C}M^*(t)$, which is equal to -2%. Indeed, if we set $\mathcal{C}M^*(t) = -3\%$, we have no solution to the optimization problem even for the base date $t_0$.

<table><caption>Table 11.25: Net-zero bond portfolio (Example 62)</caption><thead><tr><th>t</th><th>b(t<sub>0</sub>)</th><th>t<sub>0</sub></th><th>t<sub>0</sub> + 1</th><th>t<sub>0</sub> + 2</th><th>t<sub>0</sub> + 3</th><th>t<sub>0</sub> + 4</th><th>t<sub>0</sub> + 5</th><th>t<sub>0</sub> + 10</th></tr></thead><tbody><tr><td>w<sub>1</sub><sup>*</sup></td><td>20.00</td><td>4.28</td><td>13.80</td><td>20.48</td><td>26.34</td><td>19.02</td><td></td><td></td></tr><tr><td>w<sub>2</sub><sup>*</sup></td><td>19.00</td><td>34.78</td><td>38.94</td><td>42.72</td><td>46.23</td><td>49.01</td><td></td><td></td></tr><tr><td>w<sub>3</sub><sup>*</sup></td><td>17.00</td><td>21.03</td><td>13.86</td><td>7.73</td><td>2.11</td><td>0.00</td><td></td><td></td></tr><tr><td>w<sub>4</sub><sup>*</sup></td><td>13.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>No feasible solution</td><td></td></tr><tr><td>w<sub>5</sub><sup>*</sup></td><td>12.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td></td><td></td></tr><tr><td>w<sub>6</sub><sup>*</sup></td><td>8.00</td><td>39.91</td><td>33.40</td><td>29.07</td><td>25.32</td><td>31.97</td><td></td><td></td></tr><tr><td>w<sub>7</sub><sup>*</sup></td><td>6.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td></td><td></td></tr><tr><td>w<sub>8</sub><sup>*</sup></td><td>5.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td></td><td></td></tr><tr><td>AS(w)</td><td>0.00</td><td>51.72</td><td>45.34</td><td>45.27</td><td>50.89</td><td>53.98</td><td></td><td></td></tr><tr><td>MD(w)</td><td>5.48</td><td>5.48</td><td>5.48</td><td>5.48</td><td>5.48</td><td>5.48</td><td></td><td></td></tr><tr><td>DTS(w)</td><td>301.05</td><td>236.99</td><td>202.30</td><td>173.29</td><td>146.83</td><td>141.34</td><td></td><td></td></tr><tr><td>D(w | b)</td><td>0.00</td><td>0.87</td><td>0.95</td><td>1.09</td><td>1.28</td><td>1.48</td><td></td><td></td></tr><tr><td>CI(w)</td><td>160.57</td><td>112.40</td><td>104.53</td><td>97.22</td><td>90.41</td><td>84.08</td><td></td><td></td></tr><tr><td>R(w | b)</td><td>0.00</td><td>30.00</td><td>34.90</td><td>39.46</td><td>43.70</td><td>47.64</td><td></td><td></td></tr><tr><td>CM(t, w)</td><td>-1.66</td><td>-2.81</td><td>-2.57</td><td>-2.35</td><td>-2.15</td><td>-2.01</td><td></td><td></td></tr><tr><td>GI(t, w)</td><td>15.99</td><td>31.98</td><td>31.98</td><td>32.37</td><td>32.80</td><td>35.52</td><td></td><td></td></tr></tbody></table>

The reduction rate, weights, carbon momentum, green intensity and active share metrics are expressed in %, the MD metrics are measured in years, and the DTS metrics are calculated in bps.

## Empirical results

The following empirical results are taken from Barahhou et al. (2022) and Ben Slimane et al. (2023b). First, we examine the dynamic decarbonization of equity portfolios. In Figure 11.20, we show the relationship between time and tracking error volatility with respect to Scope emissions when
<figure><img src="image_514.png" /><figcaption>Figure 11.21: Tracking error volatility of dynamic decarbonized portfolios (MSCI World, June 2022, $C_3$ (0, 10, 2) constraint)</figcaption></figure>

<figure><img src="image_515.png" /></figure>

Source: MSCI (2022), Trucost (2022) & Barahhou et al. (2022).
considering the CTB and PAB decarbonization pathways and the MSCI World universe. Including Scope 3 has a significant impact on the tracking risk, especially when upstream Scope 3 emissions are considered. On average, the inclusion of Scope 3 results in a multiplication of the tracking error volatility by a factor of three. Barahhou et al. (2022) reported results considering the $C_3$ (0, 2, 1) constraint, which imposes sector neutrality. They showed that the solution may not exist even before 2030 for the PAB decarbonization pathway. In order to have acceptable solutions, the ETF industry generally uses the $C_3$ (0, 10, 2) constraint (Figure 11.21).

The previous analysis deals only with the decarbonization dimension. Barahhou et al. (2022) then introduced the transition dimension and solved the following optimization problem:

$$ w^*(t) = \arg\min \frac{1}{2} (w - b(t))^{\top} \Sigma(t) (w - b(t)) \\ \text{s.t.} \quad \begin{cases} \mathcal{CI}(t, w) \le (1 - \mathcal{R}(t_0, t)) \mathcal{CI}(t_0, b(t_0)) \\ \mathcal{CM}(t, w) \le \mathcal{CM}^*(t) \\ \mathcal{GI}(t, w) \ge (1 + \mathcal{G}) \mathcal{GI}(t, b(t)) \\ w \in C_0 \cap C_3(0, 10, 2) \end{cases} $$

where $\mathcal{CM}^*(t) = -5\%$ and $\mathcal{G} = 100\%$. In Figure 11.23, we plot the relationship between time $t$ and the tracking error volatility $\sigma(w^*(t) | b(t))$, measured in bps, when considering the PAB decarbonization pathway$^{33}$. We also report the decomposition between the decarbonization and transition dimensions. The results of these simulations clearly show that the transition dimension induces significant and additional costs. On average, we observe that the additional cost of the tracking error for the years 2022–2030 is 27, 25, 21 and 19 bps for Scopes $\mathcal{SC}_1$, $\mathcal{SC}_{1-2}$, $\mathcal{SC}_{1-3}^\text{up}$ and $\mathcal{SC}_{1-3}$, respectively. Moreover, there may be no solution to the optimization problem by 2050, especially if the carbon footprint is based on upstream/downstream Scope 3 emissions. Of course, all these results are very sensitive to the choice of the green multiplier $\mathcal{G}$ and the carbon threshold $\mathcal{CM}^*$.

The previous results are valid for the MSCI World index, which is a large investment universe with more than 1 500 stocks. Let us focus on smaller investment universes by considering the MSCI EMU and USA indexes. The results are shown in Figures 11.24 and 11.25. The tracking error volatilities for smaller universes become larger in fewer years than for the MSCI World index, and we also fail to find solutions sooner. We could separate these results by putting the Scope 1 and 2 alignment on one side and Scope 3 on the other. Looking at Scopes 1 and 2, we see that in both universes the aligned portfolio breaks earlier than the MSCI World. However, even though the MSCI EMU universe is smaller than the MSCI USA universe, we can find solutions for a longer period of time. This is due to the distribution of green revenues and carbon dynamics, which are easier to reconcile with the intensity reduction constraint for the EMU. The inclusion of Scope 3 intensities paints a different picture. Although the EMU net-zero portfolios have lower tracking errors than the USA net-zero portfolios, larger universes tend to provide longer solutions. The fact that we are unable to align the EMU portfolio in terms of Scope 3 carbon intensities after 2040 therefore highlights the difficulty of portfolio alignment for a relatively small investment universe.

In practice, many constraints can be used to construct net-zero portfolios. We have seen above that the cost of tracking error can be significant and that the solution may not exist for long time horizons. Because the net-zero portfolio excludes some assets, it may be more concentrated than the benchmark. Therefore, we may face not only diversification risk, but also liquidity risk. These risks will be reduced as the economy decarbonizes in the coming years. However, we are not immune to

<footnote>$^{33}$To assess the impact of the weight/sector constraint $C_3$ (0, 10, 2), we report the results based on the $C_0$ constraint only in Figure 11.22.</footnote>
<figure><img src="image_516.png" /><figcaption>Figure 11.23: Tracking error volatility of net-zero portfolios (MSCI World, June 2022, $C_3(0, 10, 2)$ constraint, $G = 100\%$, $CM^* = -5\%$, PAB)</figcaption></figure>

<figure><img src="image_517.png" /></figure>

Source: MSCI (2022), Trucost (2022) & Barahhou et al. (2022).
<figure><img src="image_518.png" /><figcaption>Figure 11.25: Tracking error volatility of net-zero portfolios (MSCI USA, Jun. 2022, $C_3(0, 10, 2)$ constraint, $G = 100\%$, $CM^* = -5\%$, PAB)</figcaption></figure>

<figure><img src="image_519.png" /></figure>

Source: MSCI (2022), Trucost (2022) & Barahou et al. (2022).
the possibility that carbon emissions will continue to rise in the short term. In this case, solutions will be very sensitive to the gap between the carbon target of net-zero portfolios and the carbon footprint of the economy. To illustrate the shrinkage risk of the investment universe, we calculate the number of stocks selected per sector for each optimized portfolio and divide these figures by the corresponding number of stocks in the index<fnref n="34" />. In the case of Scope **SC₁⁻₃**, the radar chart of these frequencies is shown in Figure 11.26. We observe that the investment universe is shrunk at the first date. The green area represents the removed part by 2030. With the exception of the communication services, financials, health care, information technology and real estate sectors, the investment in the other sectors is concentrated on few stocks. This shrinkage effect is also observed for small investment universes. By construction, the shrinkage of the investment universe worsens if we add other constraints. For instance, the impact of the momentum exclusion constraint is illustrated in Figure 11.27. In this case, we complete the set of constraints by the exclusion constraint {**CM**ᵢ(t) ≥ 0 ⇒ wᵢ = 0}, meaning that we exclude issuers with a positive carbon trend. We notice that the investment universe is highly reduced even from the first year. These results show that we cannot reduce the cost of net-zero investing to the cost of tracking risk. As seen above, there is also a cost of diversification risk. There is also a liquidity risk, as illustrated by Barahhou et al. (2022). Indeed, they showed that the repartition between large, mid and small caps changes. In some particular cases, they observed that the allocation to small- and micro-cap buckets increases over time.

<figure><img src="image_520.png" /><figcaption>Figure 11.26: Radar chart of investment universe shrinkage (MSCI World, June 2022, C₃ (0, 10, 2) constraint, G = 100%, CM* = -5%, PAB, Scope SC₁⁻₃)</figcaption></figure>

Source: MSCI (2022), Trucost (2022) & Barahhou et al. (2022).

<footnote n="34">³⁴For example, if the frequency for the energy sector is 25%, this means that the optimized portfolio selected 25% of the energy stocks and removed 75% of the energy investment universe.</footnote>
The results obtained by Barahhou et al. (2022) are based on simulations and do not take into account all the investment constraints that may be encountered in a net-zero portfolio, such as ESG constraints and the exclusion of some activities. Therefore, the estimated cost of tracking error risk can be viewed as a lower bound on the cost of real net-zero investing.

<figure><img src="image_521.png" /><figcaption>Figure 11.27: Impact of momentum exclusion on universe shrinkage (MSCI World, June 2022, $C_3(0, 10, 2)$ constraint, $G = 100\%$, $CM^* = -5\%$, PAB, Scope $SC_{1-3}$, $CM^+ = 0\%$)</figcaption></figure>

Source: MSCI (2022), Trucost (2022) & Barahhou et al. (2022).

We now turn to the case of corporate bonds. Using the CTB and PAB decarbonization scenarios, we get the results in Figures 11.28 and 11.29. The DTS risk is not significant and is less than 6 bps until 2030. This is not the case for the active share risk, which can reach 20% in 2030 for the PAB decarbonization pathway. When we include the transition constraints, the additional cost seems relatively low compared to what we have observed for equity investment universes. For example, Barahhou et al. (2022) found that the DTS tracking risk and active share increase by less than 1 bp and 1%, respectively, when $G$ is set to 100%. However, they found that the exclusion constraint can significantly increase costs. For example, requiring issuers to have a negative carbon momentum increases the tracking risk by 20% on average.

The case of government bonds is studied by Barahhou et al. (2023). As with corporate bonds, the integrated approach for sovereign bonds consists of several steps: (1) we need to define the decarbonization scenario at the country level; (2) we can assess the self-decarbonization of a country by considering the government's credible commitments and decarbonization plans towards a low-carbon economy; (3) a specific green intensity measure needs to measure the country's contribution to the climate transition and its greenness. Unlike corporate bonds, there are many options for choosing the net-zero scenario. We have already discussed some of them in Section 8.3.3 on page 781. We can choose between the IPCC scenarios, the NGFS scenarios, the IEA scenarios, etc. Barahhou et
<figure><img src="image_522.png" /><figcaption>Figure 11.29: Active share of dynamically decarbonized portfolios (Global Corporate, June 2022)</figcaption></figure>

<figure><img src="image_523.png" /></figure>

Source: ICE (2022), Trucost (2022) & Barahhou et al. (2022).
al. (2023) explained that the stated policies scenario (STEPS) and the announced pledges scenario (APS) cannot be used because they are not consistent with limiting global warming to 1.5°C. It is better to use the NZE scenarios from NGFS or IEA. Figure 11.30 shows the decarbonization pathway derived from the IEA NZE scenario (IEA, 2021, Figure 2.2, page 53). We distinguish between advanced economies, and emerging market and developing economies. Using 2020 as a baseline, the decarbonization pathways in terms of carbon emissions are similar for developed and developing countries. This is not the case when looking at carbon intensity, which corresponds to CO₂ emissions per capita.

<figure><img src="image_524.png" /><figcaption>Figure 11.30: IEA decarbonization pathways</figcaption></figure>

Choosing the scope of emissions and how to measure the carbon footprint is another important issue. As seen in Section 8.4.4 on page 836, we can define carbon emissions at the country level by considering a production-based or consumption-based inventory. If we prefer a carbon intensity measure, the normalization variable can be population, GDP, or public debt. Barahhou et al. (2023) showed that the distribution of carbon intensity metrics across countries is very different for these three measures. For the green footprint, they listed several metrics: government spending data<fnref n="35" /> from the IEA, government spending on environmental protection<fnref n="36" /> from the IMF's Climate Change Dashboard, the amount of green bonds issued, and environmental taxes<fnref n="37" /> by country from the IMF's

<footnote n="35">³⁵The website is https://www.iea.org/reports/government-energy-spending-tracker-2.</footnote>

<footnote n="36">³⁶This database contains the following seven time series: (1) Environmental protection expenditure, (2) Biodiversity and landscape protection expenditure, (3) Expenditure on environmental protection n.e.c., (4) Environmental protection R&D expenditure, (5) Pollution abatement expenditure, (6) Waste management expenditure and (7) Waste water management expenditure. The website for downloading the data is https://climatedata.imf.org/datasets/d22a6decd9b147fd9040f793082b219b_0/explore.</footnote>

<footnote n="37">³⁷This database contains the following five time series: (1) Environmental taxes, (2) Taxes on energy (including fuels for transport), (3) Taxes on pollution, (4) Taxes on resources and (5) Taxes on transport (excluding fuels for transport). The website for downloading the data is https://climatedata.imf.org/datasets/</footnote>
Climate Change Dashboard. Using data on carbon emissions and commitments, they also derived four types of forward-looking metrics: (1) carbon trend, (2) nationally determined contribution (NDC), (3) NDC ambition, and (4) NDC fulfillment. In particular, Barahhou et al. (2023) defined two simple criteria:

* Commitments aligned with the NZE scenario (CAS) imply that countries have NDC target emissions that fall within a certain range around the value projected by the NZE scenario;

* Emissions on track with commitments (EOTC) indicate countries with historical trends in line with their NDC commitments.

These two metrics measure the ambition and credibility of the *PAC* framework. In Table 11.26, we report the CAS and EOTC statistics found by Barahhou et al. (2023). Only 6.6% of countries have sufficient ambition and credibility consistent with the NZE scenario, while about 70% of countries do not meet these two criteria.

<table><caption>Table 11.26: CAS and EOTC statistics on Bloomberg Global Aggregate Treasuries</caption><thead><tr><th>CAS</th><th>✓</th><th>✓</th><th></th><th>No criteria</th></tr><tr><th>EOTC</th><th>✓</th><th></th><th>✓</th><th>met</th></tr></thead><tbody><tr><th>Frequency</th><td>6.6%</td><td>23.3%</td><td>0.7%</td><td>69.5%</td></tr></tbody></table>

Source: Barahhou et al. (2023, Figure 14, page 23).

Very quickly, there is no solution to the optimization problem after 2030, which means that the current ambitions of countries are not sufficient to build a net-zero sovereign bond portfolio in the long run. In Table 11.27, we report the authors' estimated first year of country exit for a given set of parameters and the GHG/GDP intensity metric. Some countries are removed from the portfolio in 2024 such as Canada, Indonesia, New Zealand and South Korea, but most of developed countries exit the portfolio in 2029. Changing the green constraint could delay the first year of exit by a year or two, but the problem is the constraint of the decarbonization pathway, which is impossible to manage after 2032.

<table><caption>Table 11.27: First year of country exit from the NZE investment portfolio (GHG/GDP intensity metric)</caption><tr><td>Australia</td><td>2025</td><td>Finland</td><td>2029</td><td>Lithuania</td><td>2025</td><td>Romania</td><td>2029</td></tr><tr><td>Austria</td><td>2029</td><td>France</td><td>2029</td><td>Luxembourg</td><td>2029</td><td>Singapore</td><td>2029</td></tr><tr><td>Belgium</td><td>2028</td><td>Germany</td><td>2029</td><td>Mexico</td><td>2029</td><td>Slovakia</td><td>2025</td></tr><tr><td>Canada</td><td>2024</td><td>Hong Kong</td><td>2029</td><td>Malaysia</td><td>2028</td><td>Slovenia</td><td>2028</td></tr><tr><td>Chile</td><td>2029</td><td>Hungary</td><td>2029</td><td>Malta</td><td>2029</td><td>South Korea</td><td>2024</td></tr><tr><td>China</td><td>2028</td><td>Indonesia</td><td>2024</td><td>Netherlands</td><td>2029</td><td>Spain</td><td>2028</td></tr><tr><td>Colombia</td><td>2029</td><td>Ireland</td><td>2029</td><td>Norway</td><td>2029</td><td>Switzerland</td><td>2029</td></tr><tr><td>Cyprus</td><td>2029</td><td>Israel</td><td>2029</td><td>New Zealand</td><td>2024</td><td>Sweden</td><td>2029</td></tr><tr><td>Czechia</td><td>2024</td><td>Italy</td><td>2029</td><td>Peru</td><td>2029</td><td>Thailand</td><td>2025</td></tr><tr><td>Denmark</td><td>2029</td><td>Japan</td><td>2029</td><td>Poland</td><td>2029</td><td>United Kingdom</td><td>2029</td></tr><tr><td>Estonia</td><td>2025</td><td>Latvia</td><td>2028</td><td>Portugal</td><td>2028</td><td>United States</td><td>2028</td></tr></table>

Source: Barahhou et al. (2023, Table 9, page 26).

<footnote>3fb1ed30d3394574b3145246846023b1_0/explore.</footnote>
<table><caption>Table 11.28: Country exclusion year by intensity metric</caption><thead><tr><th rowspan="2">Metric</th><th>GHG</th><th>GHG</th><th>CO<sub>2</sub> (production)</th><th>CO<sub>2</sub> (consumption)</th></tr><tr><th>GDP</th><th>Population</th><th>GDP</th><th>Population</th></tr></thead><tbody><tr><td>China</td><td>2028</td><td>2031</td><td>2027</td><td>2031</td></tr><tr><td>France</td><td>2029</td><td>2032</td><td>2027</td><td>2031</td></tr><tr><td>Indonesia</td><td>2024</td><td>2032</td><td>2024</td><td>2031</td></tr><tr><td>Ireland</td><td>2029</td><td>2030</td><td>2027</td><td>2030</td></tr><tr><td>Japan</td><td>2029</td><td>2032</td><td>2027</td><td>2031</td></tr><tr><td>United States</td><td>2028</td><td>2030</td><td>2026</td><td>2029</td></tr><tr><td>United Kingdom</td><td>2029</td><td>2032</td><td>2027</td><td>2031</td></tr><tr><td>Sweden</td><td>2029</td><td>2032</td><td>2027</td><td>2031</td></tr></tbody></table>

Source: Barahhou et al. (2023, Table 14, page 31).

The effect of the intensity metric is also significant. Table 11.28 shows when the country is removed from the NZE investment portfolio. We observe significant differences when we use GHG per GDP or GHG per capita. In addition, the use of a production-based or a consumption-based inventory also affects the results<fnref n="38" />. However, regardless of the metrics used, we generally found that a net-zero sovereign bond portfolio tended to overweight European countries.

The previous empirical results for equities, corporate bonds, and sovereign bonds suggest the following lessons. First, the solution is parameter and data sensitive. In particular, we need to be careful in choosing the carbon scope metric to assess the decarbonization rate. A net-zero investment policy only makes sense for a closed system. Therefore, Scope 3 emissions need to be taken into account to align a portfolio with a net-zero scenario. The problem is that we see a lack of data reliability on Scope 3 emissions today. Similarly, the solution is highly dependent on the green intensity target and the level of self-decarbonization we want to achieve. Then we have to be careful because there may be no solution to the optimization problem in the medium term. The question of no solution depends on the relative speed of the portfolio's decarbonization path relative to the economy's decarbonization path and the initial starting point. The second key finding is that portfolio alignment (or decarbonization) and net-zero construction lead to different solutions. In particular, decarbonizing a portfolio is easier than constructing a net-zero portfolio. We find that decarbonizing along CTB or PAB pathways never leads to exploding tracking errors by 2030. In fact, the real problem with decarbonization is the diversification and liquidity risk that an investor may face. These results are amplified when we add the transition dimension to the optimization program. In addition to higher tracking risk, there is no guarantee that there will always be a solution. Moreover, the introduction of the transition pillar highlights the difficulty of choosing an appropriate set of constraints for net-zero portfolios, as some metrics may be negatively correlated with others. Portfolio decarbonization is systematically a strategy that is long financial issuers and short energy, materials and utilities issuers. Therefore, we have a situation where the transition dimension of a decarbonized portfolio is weaker than that of the benchmark portfolio, as green solutions are also located in carbon-intensive sectors. It is therefore crucial to distinguish between issuers with a high carbon footprint that will not participate in the transition and those that will reduce their carbon emissions and find low-carbon solutions. The third key finding is that portfolio

<footnote n="38">³⁸In this case, estimated metrics are generally available for CO₂ emissions than for all GHG emissions because consumption-based estimates are obtained with input-output analysis. In addition, it is traditionally accepted that the intensity measure uses GDP as the normalization variable for the production-based inventory, because it depends on national production, and the population for the consumption-based inventory, because it depends on population size.</footnote>
decarbonization and alignment are two processes of exclusion. This means that it is quite impossible to achieve net-zero alignment without allowing the algorithm to exclude companies (or countries) from the benchmark. For example, the optimization program will generally not find a solution if it imposes non-zero lower bounds. As a result, some key players in the transition, such as energy and utility companies, unfortunately disappear. Furthermore, imposing sector neutrality can lead to similar problems in finding a solution. The final lesson is that it is easier to implement net-zero in bonds than in equities. At first glance, this result may seem surprising, since there is no reason why net-zero should affect the equity and bond markets differently. In fact, there are two possible explanations. First, the structure of equity and bond indexes is different, with the latter having a more balanced allocation across sectors and a high exposure to financial issuers. Second, bond indexes are strongly influenced by new fresh capital, while equity indexes are sticky to the stock of existing capital. This is because the primary bond market is very active, which implies a significant impact on the secondary market. Indeed, bonds mature and are replaced by new, greener bonds. The primary market then helps to achieve net-zero in bonds. This is not the case in the stock market, where IPOs and capital increases are only a small part of the secondary market. This means that portfolio holdings change faster for bond indexes than for equity indexes. Therefore, the greenness of bond indexes increases faster than the greenness of equity indexes. All of these factors suggests that the cost of implementing net-zero investments relative to traditional investments will be higher for equity portfolios than for bond portfolios, and that the bond market will benefit more quickly from the transition to a low-carbon economy.

### 11.3.2 Core-satellite approach

We have seen that the comprehensive integrated approach can sometimes be difficult to implement because today, on average, carbon intensities are positively correlated with green intensities. This means that the greenness of the economy is not necessarily found in companies with low carbon footprints. Therefore, a second approach has emerged that is easier to implement. It consists of adopting a core-satellite strategy, where decarbonization is applied to the core portfolio, while the objective of the satellite portfolio is to finance the transition to a low-carbon economy. In the financial literature, the core portfolio is called the net-zero decarbonization portfolio, while the satellite portfolio is called the net-zero contribution portfolio, but other terms are used, as shown in Table 11.29.

<table><caption>Table 11.29: The two building block approach</caption><thead><tr><th>Decarbonizing the portfolio</th><th>Financing the transition</th></tr></thead><tbody><tr><td><ul><li>Net-zero decarbonization portfolio</li><li>Net-zero transition portfolio</li><li>Dynamic low-carbon portfolio</li></ul></td><td><ul><li>Net-zero contribution portfolio</li><li>Net-zero funding portfolio</li><li>Net-zero transformation portfolio</li></ul></td></tr></tbody></table>

This is equivalent to splitting the problem into two sub-problems. The goal of the first sub-problem is to decarbonize and manage the carbon footprint of the investment. The goal of the second sub-problem is to contribute to increasing the green footprint of the economy. The two sub-problems are summarized in Table 11.30. The core portfolio is more of a top-down allocation process and exclusion strategy, where the central climate risk metric is carbon intensity. The satellite
portfolio is more of a bottom-up allocation process and asset selection strategy, where the central climate risk metric is green intensity. This approach also has the advantage of making the allocation between the two net-zero strategies clear. Of course, the allocation $\alpha(t)$ to the satellite can be dynamic and change over time as the world and economy progresses towards net-zero.

<figure><img src="image_525.png" /><figcaption>Table 11.30: The core-satellite approach</figcaption></figure>

## Core portfolio

The core portfolio is the largest part of the net-zero investment portfolio and the objective is to manage the decarbonization dimension. Therefore, only carbon metrics are used to assess the carbon footprint of the portfolio and its trajectory. If we write the alignment process in terms of portfolio optimization, a typical program for the equity bucket looks like this:

$$ w^*(t) = \arg \min \frac{1}{2} (w - b(t))^\top \Sigma(t) (w - b(t)) \quad (11.24) \\ \text{s.t.} \quad \begin{cases} \mathcal{CI}(t, w) \le (1 - \mathcal{R}(t_0, t)) \mathcal{CI}(t_0, b(t_0)) \\ \mathcal{CM}(t, w) \le \mathcal{CM}^*(t) \\ 0 \le w_i \le 1, \{\mathcal{CM}_i(t) \le \mathcal{CM}^+\} \\ w \in \Omega_0 \cap \Omega \end{cases} $$

This is exactly the optimization problem defined in Equation (11.22), except that the green intensity constraints are removed because the green footprint is managed through the satellite portfolio. For the bond bucket, we get a similar optimization problem:

$$ w^*(t) = \arg \min_{\mathcal{D}} (w | b(t)) \quad (11.25) \\ \text{s.t.} \quad \begin{cases} \mathcal{CI}(t, w) \le (1 - \mathcal{R}(t_0, t)) \mathcal{CI}(t_0, b(t_0)) \\ \mathcal{CM}(t, w) \le \mathcal{CM}^*(t) \\ 0 \le w_i \le 1, \{\mathcal{CM}_i(t) \le \mathcal{CM}^+\} \\ w \in \Omega_0 \cap \Omega \end{cases} $$

In both cases, the implementation of the portfolio alignment process considers a global decarbonization scenario constraint, a self-decarbonization constraint, and an exclusion constraint. Ben Slimane et al. (2023b) also suggested that a number of specific sectoral decarbonization pathways can be considered, consistent with the International Energy Agency's sequencing principles. In particular, the issue of electricity is central to the net zero scenario. Electricity has to be green by 2035. From a financial perspective, this means that we will have to finance the transition of the power sector on
a massive scale. However, if we only implement a global decarbonization scenario, Ben Slimane et al. (2023b) noted that utility issuers are quickly excluded from the net-zero optimization program. The core net-zero portfolio invests in pure players that produce green power and doesn’t give others time to transform their business models. This is not consistent with the goal of net-zero investing. Ben Slimane et al. (2023b) then proposed to constrain the optimization problem to follow the NZE scenario for the electricity sector. Using the results on page 1044, the constraint to meet a reduction rate for a given sector $Sector_j$ can be expressed as:

$$ \frac{\sum_{i=1}^{n} 1 \{i \in Sector_j\} w_i CI_i}{\sum_{i=1}^{n} 1 \{i \in Sector_j\} w_i} = CI(Sector_j, R_j) $$

where $CI(Sector_j, R_j)$ is the carbon intensity target for the given sector:

$$ CI(Sector_j, R_j) = (1 - R_j) \frac{\sum_{i=1}^{n} 1 \{i \in Sector_j\} b_i CI_i}{\sum_{i=1}^{n} 1 \{i \in Sector_j\} b_i} $$

We deduce that:

$$ \sum_{i=1}^{n} 1 \{i \in Sector_j\} w_i CI_i = CI(Sector_j, R_j) \sum_{i=1}^{n} 1 \{i \in Sector_j\} w_i $$

which is equivalent to the following constraint:

$$ \sum_{i=1}^{n} 1 \{i \in Sector_j\} w_i (CI_i - CI(Sector_j, R_j)) = 0 \Leftrightarrow (s_j \circ (CI_i - CI_j^*))^\top w = 0 $$

where $CI_j^* = CI(Sector_j, R_j)$. Since this is a linear equation, the previous optimization problems with this constraint remain a QP or LP problem and can be solved easily. In Figures 11.31 and 11.32, we reproduce the tracking error volatility $\sigma(w | b)$ of equity portfolios and the active risk $D(w | b)$ of bond portfolios when we impose the IEA NZE scenario for the electricity sector, $CM^* = -3.5\%$ and $CM^+ = 10\%$, and we consider different NZE decarbonization scenarios (IEA, NZAOA, CTB and PAB). In these simulations, Ben Slimane et al. (2023b) found that if we allocate the tracking risk between the different constraints, the two most important contributors are the global decarbonisation pathway and the exclusion constraints, while the costs of self-decarbonisation and the specific electricity decarbonisation pathway are relatively low.

## Satellite portfolio

While the core portfolio aims to implement decarbonization policies, the satellite portfolio aims to finance the transition to a low-carbon economy and monitor green intensity. As the core-satellite approach is implemented in strategic asset allocation or multi-asset portfolios, the investment universe is diversified and typically consists of the following asset classes:

* Green, sustainability and sustainability-linked bonds

* Green stocks

* Green infrastructure

* Sustainable real estate
<figure><img src="image_526.png" /><figcaption>Figure 11.32: Active risk of decarbonized portfolios (Global Corporate, December 2021, $\mathcal{CM}^* = -3.5\%$, $\mathcal{CM}^+ = 10\%$, IEA NZE electricity sector scenario)</figcaption></figure>

<figure><img src="image_527.png" /></figure>

Source: Ben Slimane et al. (2023b).
<table><caption>Table 11.31: Main sub-industries of the net-zero satellite portfolio (GICS level 4)</caption><thead><tr><th>Code</th><th>Level 4</th><th>Level 2</th><th>Level 1</th></tr></thead><tbody><tr><td>15102010</td><td>Construction Materials</td><td></td><td></td></tr><tr><td>15104010</td><td>Aluminium</td><td></td><td></td></tr><tr><td>15104020</td><td>Diversified Metals & Mining</td><td></td><td></td></tr><tr><td>15104025</td><td>Copper</td><td>Materials</td><td>Materials</td></tr><tr><td>15104040</td><td>Precious Metals & Minerals</td><td></td><td></td></tr><tr><td>15104045</td><td>Silver</td><td></td><td></td></tr><tr><td>15104050</td><td>Steel</td><td></td><td></td></tr><tr><td>20102010</td><td>Building Products</td><td></td><td></td></tr><tr><td>20103010</td><td>Construction & Engineering</td><td></td><td></td></tr><tr><td>20104010</td><td>Electrical Components & Equipment</td><td>Capital Goods</td><td></td></tr><tr><td>20104020</td><td>Heavy Electrical Equipment</td><td></td><td></td></tr><tr><td>20106010</td><td>Construction Machinery & Heavy Transportation Eqpt.</td><td></td><td></td></tr><tr><td>20106015</td><td>Agricultural & Farm Machinery</td><td></td><td></td></tr><tr><td>20201050</td><td>Environmental & Facilities Services</td><td>Commercial & Professional Services</td><td>Industrials</td></tr><tr><td>20304010</td><td>Rail Transportation</td><td></td><td></td></tr><tr><td>20305010</td><td>Airport Services</td><td></td><td></td></tr><tr><td>20305020</td><td>Highways & Railtracks</td><td>Transportation</td><td></td></tr><tr><td>20305030</td><td>Marine Ports & Services</td><td></td><td></td></tr><tr><td>25101010</td><td>Automotive Parts & Equipment</td><td>Automobiles & Components</td><td></td></tr><tr><td>25102010</td><td>Automobile Manufacturers</td><td></td><td></td></tr><tr><td>25201010</td><td>Consumer Electronics</td><td>Consumer Durables & Apparel</td><td>Consumer Discretionary</td></tr><tr><td>25201030</td><td>Homebuilding</td><td></td><td></td></tr><tr><td>25201040</td><td>Household Appliances</td><td></td><td></td></tr><tr><td>30202010</td><td>Agricultural Products</td><td>Services Food, Beverage & Tobacco</td><td>Consumer Staples</td></tr><tr><td>55101010</td><td>Electric Utilities</td><td></td><td></td></tr><tr><td>55103010</td><td>Multi-Utilities</td><td></td><td></td></tr><tr><td>55104010</td><td>Water Utilities</td><td>Utilities</td><td>Utilities</td></tr><tr><td>55105020</td><td>Renewable Electricity</td><td></td><td></td></tr><tr><td>60201030</td><td>Real Estate Development</td><td>Real Estate Management & Development</td><td>Real Estate</td></tr></tbody></table>

Source: Ben Slimane et al. (2023b, Table 14, page 44).
The net-zero transition analysis on page 1031 shows that the list of major sectors to be financed is relatively small. The mapping to level 4 of the GICS classification is shown in Table 11.31. In fact, we identify 29 sub-industries out of the 163 included in the GICS classification. Looking at the major sub-industries that can be included in the net-zero satellite portfolio, four of the eleven GICS sectors are over-represented, two are included, and the other five are excluded. The four over-represented sectors are Industrials, Materials, Consumer Discretionary and Utilities. The Industrials sector is divided into three industry groups: Capital Goods, which includes sub-industries related to machinery, equipment and construction; Transportation, where we look at how to improve public transportation systems and their infrastructure; and finally, Commercial & Professional Services, which includes the Environmental & Facilities Services sub-industry, which mainly covers waste management and pollution control. The Materials sector is characterized by the various materials used in the energy transition, such as aluminium, copper, steel, etc. The Consumer Discretionary sector is divided into two industry groups, Automobiles & Components, where we find auto manufacturers and auto parts, and Consumer Durables & Apparel, related to housing and appliances. Finally, the Utilities sectors will show different types of utilities needed for the transition, such as electric utilities, water utilities, or even renewable energy. Then we include two sub-industries that belong to the Consumer Staples and Real Estate sectors, respectively: agricultural products and real estate development. In summary, this means that not all sectors are represented. Figure 11.33 shows the four levels of the GICS classification<fnref n="39" /> and indicates which sub-industry falls within the definition of the satellite investment universe.

**Remark 146** It is clear that the GICS classification is not relevant when considering a net-zero investing framework. For example, there is no sector such as electricity storage, hydrogen storage, photovoltaic electricity generation, wind electricity generation, nuclear electricity generation in existing plants, etc. The NACE classification is more appropriate, and has the advantage to be in line with the EU green taxonomy. Nevertheless, GICS is the classification used by investors.

**Green bonds** We look at the Bloomberg database of the GSS+ investment universe. For each bond, Bloomberg indicates whether it is a green, social, sustainability, sustainability-linked or conventional bond. Issue amounts are shown in Table 11.32. In 2022, 1 784 green bonds were issued for a total of $531.6 bn. This represents 15% of the net-zero financing needs<fnref n="40" />. For the other categories, the amount to be issued in 2022 is equal to $152.8 bn for social bonds, $174.8 bn for sustainability bonds and $144.3 bn for sustainability-linked bonds. As explained by Ben Slimane et al. (2023a), social bonds are not net-zero transition instruments, but more conventional bonds to finance social debt and social infrastructure. Therefore, if we look at a broad definition of the net-zero fixed income universe (green, sustainability and sustainability-linked bonds), we get a total of $850.7 bn, which can be seen as the upper bound of current investment opportunities. This is less than 25% of the $3.5 tn previously required to achieve net-zero.

In Figure 11.34 we compare the performance of the Bloomberg Global Green Bond index with the performance of the Bloomberg Global Aggregate index. We see that there is a high tracking risk. Between January 2014 and September 2023, the volatility of the tracking error is 2.7%. Several factors explain this high figure: sector allocation, duration, credit risk, etc. Investors must therefore accept a higher active risk for the satellite portfolio than for the core portfolio.

<footnote n="39">³⁹It is described in Box 11.1 on page 1066.</footnote>

<footnote n="40">⁴⁰On page 10.4.1 we saw that the net-zero funding requirement is equivalent to $3.5 tn per year.</footnote>
<table><caption>Figure 11.33: Narrow specification of the satellite investment universe</caption><thead><tr><th>Sector</th><th>Industry Group</th><th>Industry</th><th>Sub-industry</th><th>Satellite</th></tr></thead><tbody><tr><td>10</td><td></td><td></td><td></td><td></td></tr><tr><td>15</td><td></td><td></td><td></td><td></td></tr><tr><td>20</td><td></td><td></td><td></td><td></td></tr><tr><td>25</td><td></td><td></td><td></td><td></td></tr><tr><td>30</td><td></td><td></td><td></td><td></td></tr><tr><td>35</td><td></td><td></td><td></td><td></td></tr><tr><td>40</td><td></td><td></td><td></td><td></td></tr><tr><td>45</td><td></td><td></td><td></td><td></td></tr><tr><td>50</td><td></td><td></td><td></td><td></td></tr><tr><td>55</td><td></td><td></td><td></td><td></td></tr><tr><td>60</td><td></td><td></td><td></td><td></td></tr></tbody></table>

Source: Ben Slimane et al. (2023b, Figure 31, page 45).

<table><caption>Table 11.32: GSS+ bond issuance</caption><thead><tr><th rowspan="2">Year</th><th colspan="2">Green</th><th colspan="2">Social</th><th colspan="2">Sustainability</th><th colspan="2">SLB</th></tr><tr><th>#</th><th>$ bn</th><th>#</th><th>$ bn</th><th>#</th><th>$ bn</th><th>#</th><th>$ bn</th></tr></thead><tbody><tr><td>2022</td><td>1 784</td><td>531.6</td><td>542</td><td>152.8</td><td>614</td><td>174.8</td><td>382</td><td>144.3</td></tr><tr><td>2021</td><td>1 971</td><td>686.1</td><td>554</td><td>242.1</td><td>646</td><td>233.2</td><td>343</td><td>161.5</td></tr><tr><td>2020</td><td>1 076</td><td>291.2</td><td>273</td><td>172.0</td><td>308</td><td>154.8</td><td>47</td><td>16.5</td></tr><tr><td>2019</td><td>877</td><td>268.0</td><td>99</td><td>22.2</td><td>333</td><td>85.2</td><td>18</td><td>8.9</td></tr><tr><td>2018</td><td>582</td><td>165.3</td><td>48</td><td>16.5</td><td>52</td><td>22.1</td><td>1</td><td>2.2</td></tr><tr><td>2017</td><td>472</td><td>160.9</td><td>46</td><td>11.8</td><td>17</td><td>9.2</td><td>1</td><td>0.2</td></tr><tr><td>2016</td><td>285</td><td>99.7</td><td>14</td><td>2.2</td><td>16</td><td>6.6</td><td>0</td><td>0.0</td></tr></tbody></table>

Source: Bloomberg (2023), GSS+ Instrument Indicator & Author's calculations.
<figure><img src="image_528.png" /><figcaption>Figure 11.34: Performance and duration of the Bloomberg Global Green Bond and Aggregate indexes</figcaption></figure>

Source: Bloomberg (2023) & Author's calculations.

**Green stocks** We can invest in green stocks through thematic funds, which focus on a specific theme, or through an equity basket. In general, the second approach is used by equity managers, while the first is preferred by other fund managers, especially multi-asset managers. The emergence of economic, social and technological megatrends is strongly influencing the expansion of thematic funds. Some funds address the environmental challenge with a broader range of climate policies and solutions, investing only in companies with the best environmental practices and targets. These funds do not necessarily focus on net-zero transition. However, they generally include renewable energy stocks. Some more specific thematic funds have been developed to address net-zero issues such as clean energy, hydrogen, water management and future mobility. Figure 11.35 shows the performance of four thematic equity indexes: Bloomberg BioEnergy, Bloomberg Hydrogen, MSCI Future Mobility and MSCI New Energy. These indexes have a high risk of tracking error relative to the MSCI World index. On average, the tracking error volatility of the satellite equity portfolio is expected to be around 20%. The fund manager can also develop a stock-picking process related to the net-zero theme. In this case, he can build a screening based on green revenue share, green capex or green opex measures. More generally, the idea behind building a satellite equity portfolio with a basket of stocks is to select stocks according to a green intensity measure that can be aligned with the EU taxonomy.

**Green infrastructure** The European Commission<fnref n="41" /> defines green infrastructure as “a strategically planned network of natural and semi-natural areas with other environmental features, designed and managed to deliver a wide range of ecosystem services, while also enhancing biodiversity.” Green

<footnote n="41">⁴¹https://environment.ec.europa.eu/topics/nature-and-biodiversity/green-infrastructure_en.</footnote>
<figure><img src="image_529.png" /><figcaption>Figure 11.35: Performance and tracking error volatility of thematic equity indexes</figcaption></figure>

Source: Bloomberg (2023), MSCI (2023) & Author's calculations.

infrastructure is implemented in a variety of sectors, from energy through energy transmission in-
frastructure, water through natural water retention measures or sustainable urban drainage systems,
to the urban landscape with street trees to help sequester carbon or green roofs to help regulate
the temperature of buildings. The cost of implementing green infrastructure is in the identification,
mapping, planning and creation of the infrastructure, but the environmental, economic and social
benefits make it worthwhile. Funds that assess infrastructure needs are emerging in the market and
typically invest in owners of sustainable infrastructure assets as well as companies that are leaders in
infrastructure investment. In addition to infrastructure funds, investors are also considering direct
investments such as green car parks, water infrastructure and flood defences.

**Sustainable real estate** The real estate sector emits significant amounts of CO2 through building operations, building materials, and construction. Action is needed in the construction of new buildings, but also in the renovation of existing buildings. For existing buildings, it is very important to reduce energy consumption, eliminate emissions from energy and refrigerants, and reduce or eliminate the use of fossil fuels. This is done by improving equipment such as insulation, ventilation, and the use of renewable energy, as well as optimizing operations by installing GHG monitors or adjusting temperature settings. New construction must be energy and carbon efficient, taking into account new and clean technologies. Sustainable real estate funds have entered the market, typically targeting multiple sectors and countries with a specific allocation to achieve net-zero by 2050. They mostly follow the CRREM (carbon risk real estate monitor) pathway, targeting 1.5°C/2°C using a Paris-aligned decarbonization pathway per country and building type, ranging from office buildings to retail stores and hotels.
Allocation process

In the following, we assume that the core-satellite portfolio is invested in stocks and bonds, but the inclusion of green infrastructure and sustainable real estate is straightforward.

**The stock/bond mix allocation** Let $\alpha_{equity}$ and $\alpha_{bond}$ be the proportions of stocks and bonds in the multi-asset portfolio. Let $\alpha^{satellite}$ be the weight of the satellite portfolio. The core allocation is given by the vector $(\alpha_{equity}^{core}, \alpha_{bond}^{core})$, while the satellite allocation is defined by $(\alpha_{equity}^{satellite}, \alpha_{bond}^{satellite})$. We have the following identities:

$$ \left\{ \begin{aligned} \alpha_{equity} &= (1 - \alpha^{satellite}) \alpha_{equity}^{core} + \alpha^{satellite} \alpha_{equity}^{satellite} \\ \alpha_{bond} &= (1 - \alpha^{satellite}) \alpha_{bond}^{core} + (1 - \alpha^{satellite}) \alpha_{bond}^{satellite} \end{aligned} \right. \qquad (11.26) $$

In general, the fund manager targets a strategic asset allocation at the portfolio level, i.e., the proportions $\alpha_{equity}$ and $\alpha_{bond}$ are given. For example, a defensive portfolio corresponds to a 20/80 constant mix strategy, while the 50/50 allocation is known as a balanced portfolio. Another famous allocation rule is the 60/40 portfolio, which is 60% in stocks and 40% in bonds. A first solution to Equation (11.26) is to maintain the same proportion of stocks and bonds in the core and satellite portfolios:

$$ \left\{ \begin{aligned} \alpha_{equity}^{core} &= \alpha_{equity}^{satellite} = \alpha_{equity} \\ \alpha_{bond}^{core} &= \alpha_{bond}^{satellite} = \alpha_{bond} \end{aligned} \right. $$

However, this solution is not always satisfactory. In fact, the satellite portfolio generally has more bonds than the overall portfolio because the investment universe of green bonds is larger than the investment universe of green stocks. As a result, bonds are overweighted in the satellite portfolio. A second solution to Equation (11.26) is to calculate the proportion of bonds in the core portfolio relative to the proportion of bonds in the satellite portfolio:

$$ \alpha_{bond}^{core} = \frac{\alpha_{bond} - \alpha^{satellite} \alpha_{bond}^{satellite}}{1 - \alpha^{satellite}} $$

**Example 63** We consider a 60/40 constant mix strategy. The satellite portfolio represents 10% of the net zero investments. We assume that the satellite portfolio has 70% exposure to green bonds.

We have $\alpha_{equity} = 60\%$, $\alpha_{bond} = 40\%$, $\alpha^{core} = 90\%$, $\alpha^{satellite} = 10\%$ and $\alpha_{bond}^{satellite} = 70\%$. We deduce that:

$$ \alpha_{bond}^{core} = \frac{0.40 - 0.10 \times 0.70}{1 - 0.10} = \frac{33}{90} = 36.67\% $$

The core allocation is then (63.33%, 36.67%), while the satellite allocation is (30%, 70%). We check that:

$$ \left\{ \begin{aligned} \alpha_{equity} &= 0.90 \times \left(1 - \frac{33}{90}\right) + 0.10 \times 0.30 = 60\% \\ \alpha_{bond} &= 0.90 \times \frac{33}{90} + 0.10 \times 0.70 = 40\% \end{aligned} \right. $$

In Table 11.33, we report the values taken by $\alpha_{bond}^{core}$ when we target several constant mix strategies (60/40, 50/50, and 20/80). We assume that $\alpha_{bond}^{satellite}$ is set to 70% or 90% and consider several values of $\alpha^{satellite}$. For example, if green bonds account for 70% of the satellite portfolio's allocation and the satellite has a weight of 25%, the core portfolio's bond allocation must be set to 30% to achieve a 60/40 constant mix strategy.
<table><caption>Table 11.33: Calculating the bond allocation in the core portfolio ($\alpha_{\text{bond}}^{\text{core}}$ in %)</caption><thead><tr><th>Strategy</th><th colspan="3">60/40</th><th colspan="3">50/50</th><th colspan="3">20/80</th></tr><tr><th>α<sup>satellite</sup><sub>bond</sub></th><th>70.0</th><th>80.0</th><th>90.0</th><th>70.0</th><th>80.0</th><th>90.0</th><th>70.0</th><th>80.0</th><th>90.0</th></tr></thead><tbody><tr><th>0%</th><td>40.0</td><td>40.0</td><td>40.0</td><td>50.0</td><td>50.0</td><td>50.0</td><td>80.0</td><td>80.0</td><td>80.0</td></tr><tr><th>1%</th><td>39.7</td><td>39.6</td><td>39.5</td><td>49.8</td><td>49.7</td><td>49.6</td><td>80.1</td><td>80.0</td><td>79.9</td></tr><tr><th>5%</th><td>38.4</td><td>37.9</td><td>37.4</td><td>48.9</td><td>48.4</td><td>47.9</td><td>80.5</td><td>80.0</td><td>79.5</td></tr><tr><th>α<sup>satellite</sup></th><td>10%</td><td>36.7</td><td>35.6</td><td>34.4</td><td>47.8</td><td>46.7</td><td>45.6</td><td>81.1</td><td>80.0</td><td>78.9</td></tr><tr><th>15%</th><td>34.7</td><td>32.9</td><td>31.2</td><td>46.5</td><td>44.7</td><td>42.9</td><td>81.8</td><td>80.0</td><td>78.2</td></tr><tr><th>20%</th><td>32.5</td><td>30.0</td><td>27.5</td><td>45.0</td><td>42.5</td><td>40.0</td><td>82.5</td><td>80.0</td><td>77.5</td></tr><tr><th>25%</th><td>30.0</td><td>26.7</td><td>23.3</td><td>43.3</td><td>40.0</td><td>36.7</td><td>83.3</td><td>80.0</td><td>76.7</td></tr></tbody></table>

**Tracking error risk of the core-satellite portfolio** Let $w$, $w^{\text{core}}$, $w^{\text{satellite}}$, and $b$ be the core-satellite, core, satellite and benchmark portfolios, respectively. The return of the core-satellite portfolio is equal to:

$$R(w) = (1 - \alpha^{\text{satellite}}) \underbrace{(\alpha_{\text{equity}}^{\text{core}} R(w_{\text{equity}}^{\text{core}}) + \alpha_{\text{bond}}^{\text{core}} R(w_{\text{bond}}^{\text{core}}))}_{\text{Core portfolio's return}} + \\ \underbrace{\alpha^{\text{satellite}} (\alpha_{\text{equity}}^{\text{satellite}} R(w_{\text{equity}}^{\text{satellite}}) + \alpha_{\text{bond}}^{\text{satellite}} R(w_{\text{bond}}^{\text{satellite}}))}_{\text{Satellite portfolio's return}}$$

By construction, we have $\alpha_{\text{equity}}^{\text{core}} + \alpha_{\text{bond}}^{\text{core}} = 1$ and $\alpha_{\text{equity}}^{\text{satellite}} + \alpha_{\text{bond}}^{\text{satellite}} = 1$. The proportion invested in equities is equal to $\alpha_{\text{equity}} = (1 - \alpha^{\text{satellite}}) \alpha_{\text{equity}}^{\text{core}} + \alpha_{\text{equity}}^{\text{satellite}} \alpha_{\text{bond}}^{\text{satellite}}$ whereas the proportion invested in bonds is the complementary part ($\alpha_{\text{bond}} = 1 - \alpha_{\text{equity}}$). We deduce that:

$$R(w) = \tilde{\alpha}^{\top} \tilde{R}(w)$$

where $\tilde{R}(w) = (R(w_{\text{equity}}^{\text{core}}), R(w_{\text{bond}}^{\text{core}}), R(w_{\text{equity}}^{\text{satellite}}), R(w_{\text{bond}}^{\text{satellite}}))$ and:

$$\tilde{\alpha} = \begin{pmatrix} (1 - \alpha^{\text{satellite}}) \alpha_{\text{equity}}^{\text{core}} \\ (1 - \alpha^{\text{satellite}}) \alpha_{\text{bond}}^{\text{core}} \\ \alpha_{\text{equity}}^{\text{satellite}} \alpha_{\text{equity}}^{\text{satellite}} \\ \alpha_{\text{bond}}^{\text{satellite}} \alpha_{\text{bond}}^{\text{satellite}} \end{pmatrix}$$

The benchmark portfolio's return is given by:

$$R(b) = \alpha_{\text{equity}} R(b_{\text{equity}}) + (1 - \alpha_{\text{equity}}) R(b_{\text{bond}})$$

Ben Slimane et al. (2023b) showed that:

$$R(b) = (1 - \alpha^{\text{satellite}}) (\alpha_{\text{equity}}^{\text{core}} R(b_{\text{equity}}) + \alpha_{\text{bond}}^{\text{core}} R(b_{\text{bond}})) + \\ \alpha^{\text{satellite}} (\alpha_{\text{equity}}^{\text{satellite}} R(b_{\text{equity}}) + \alpha_{\text{bond}}^{\text{satellite}} R(b_{\text{bond}}))$$

We deduce that:

$$R(b) = \tilde{\alpha}^{\top} \tilde{R}(b)$$

where $\tilde{R}(b) = (R(b_{\text{Equity}}), R(b_{\text{bond}}), R(b_{\text{equity}}), R(b_{\text{bond}}))$. The tracking error is defined as:

$$e = R(w) - R(b) = \tilde{\alpha}^{\top} (\tilde{R}(w) - \tilde{R}(b))$$
Let $\tilde{\Sigma}(w | b)$ be the 4 × 4 covariance matrix of $\tilde{R}(w) - \tilde{R}(b)$. We conclude that the tracking error volatility of the core-satellite portfolio has the following expression:

$$ \sigma(w | b) = \sqrt{\tilde{\alpha}^{\top} \tilde{\Sigma}(w | b) \tilde{\alpha}} = \sqrt{(\tilde{\alpha} \circ \tilde{\sigma}(w | b))^{\top} \tilde{\rho}(w | b) (\tilde{\alpha} \circ \tilde{\sigma}(w | b))} $$

where $\tilde{\rho}(w | b)$ is the correlation matrix of $R(w) - R(b)$ and $\tilde{\sigma}(w | b)$ is the vector of tracking error volatilities:

$$ \tilde{\sigma}(w | b) = \begin{pmatrix} \sigma(w_{\text{equity}}^{\text{core}} | b_{\text{equity}}) \\ \sigma(w_{\text{bond}}^{\text{core}} | b_{\text{bond}}) \\ \sigma(w_{\text{equity}}^{\text{satellite}} | b_{\text{equity}}) \\ \sigma(w_{\text{bond}}^{\text{satellite}} | b_{\text{bond}}) \end{pmatrix} $$

**Example 64** The tracking error volatilities are 2% for the core equity portfolio, 25 bps for the core bond portfolio, 20% for the satellite equity portfolio, and 3% for the satellite bond portfolio. To define the correlation matrix $\tilde{\rho}(w | b)$, we assume an 80% correlation between the two equity baskets, a 50% correlation between the two bond baskets, and a 0% correlation between the equity and bond baskets. We consider a 60/40 constant mix strategy. The satellite portfolio represents 10% of the net zero portfolio and has 70% exposure to green bonds.

We compute the tracking error covariance matrix $\tilde{\Sigma}(w | b)$ as follows: the tracking error variance for the core equity portfolio is $\tilde{\Sigma}_{1,1}(w | b) = 0.02^2$, the tracking error variance for the satellite equity portfolio is $\tilde{\Sigma}_{3,3}(w | b) = 0.20^2$, the tracking error covariance for the two core portfolios is $\tilde{\Sigma}_{1,2}(w | b) = 0 \times 0.02 \times 0.0025$, the tracking error covariance for the core equity portfolio and the satellite equity portfolio is $\tilde{\Sigma}_{1,3}(w | b) = 0.80 \times 0.02 \times 0.20$, and so on. Finally, we get:

$$ \tilde{\Sigma}(w | b) = \begin{pmatrix} 4 & 0 & 32 & 0 \\ 0 & 0.0625 & 0 & 0.375 \\ 32 & 0 & 400 & 0 \\ 0 & 0.375 & 0 & 9 \end{pmatrix} \times 10^{-4} $$

Since we have $\tilde{\alpha} = (57\%, 33\%, 3\%, 7\%)$, we deduce that $\sigma(w | b) = 1.68\%$.

We consider several constant mix strategies with the same allocation within the core and satellite portfolios: $\alpha_{\text{equity}}^{\text{core}} = \alpha_{\text{equity}}^{\text{satellite}} = \alpha_{\text{equity}}$ and $\alpha_{\text{bond}}^{\text{core}} = \alpha_{\text{bond}}^{\text{satellite}} = \alpha_{\text{bond}}$. We assume that the tracking error volatilities are 2% for the core equity portfolio, 25 bps for the core bond portfolio, 20% for the satellite equity portfolio, and 3% for the satellite bond portfolio. We consider a lower and an upper bound on the correlation matrix $\tilde{\rho}(w | b)$. For the lower bound, we set all cross-correlations to zero. For the upper bound, we assume an 80% correlation between the two equity baskets, an 80% correlation between the two bond baskets, and a 0% correlation between the equity and bond baskets. The results$^{42}$ are shown in Table 11.34 for three different values of $\alpha^{\text{satellite}}$. When the satellite weight is set to 20%, the tracking error volatility for the pure bond allocation is between 63 and 80 bps. For the 60/40 constant mix allocation, we obtain $\sigma(w | b) \in [2.60\%, 3.68\%]$, while for the pure equity allocation we have $\sigma(w | b) \in [4.31\%, 5.60\%]$. Ben Slimane et al. (2023b) estimated that the tracking error volatility of a net-zero core-satellite portfolio is currently around 3% for a 60/40 constant mix strategy.

<footnote>$^{42}$Defensive, balanced and dynamic portfolios correspond to 20/80, 50/50 and 80/20 constant mix strategies, respectively.</footnote>
<table><caption>Table 11.34: Estimation of the tracking error volatility of the core-satellite portfolio (in %)</caption><thead><tr><th></th><th>α<sup>satellite</sup></th><th>Bond</th><th>Defensive</th><th>Balanced</th><th>60/40</th><th>Dynamic</th><th>Equity</th></tr></thead><tbody><tr><td rowspan="3">Lower bound</td><td>10%</td><td>0.38</td><td>0.62</td><td>1.36</td><td>1.62</td><td>2.15</td><td>2.69</td></tr><tr><td>20%</td><td>0.63</td><td>1.00</td><td>2.18</td><td>2.60</td><td>3.45</td><td>4.31</td></tr><tr><td>30%</td><td>0.92</td><td>1.43</td><td>3.11</td><td>3.71</td><td>4.93</td><td>6.16</td></tr><tr><td rowspan="3">Upper bound</td><td>10%</td><td>0.53</td><td>1.18</td><td>2.16</td><td>2.49</td><td>3.15</td><td>3.80</td></tr><tr><td>20%</td><td>0.80</td><td>1.76</td><td>3.20</td><td>3.68</td><td>4.64</td><td>5.60</td></tr><tr><td>30%</td><td>1.07</td><td>2.34</td><td>4.24</td><td>4.87</td><td>6.13</td><td>7.40</td></tr></tbody></table>

**Dynamic allocation** The allocation between the core and satellite portfolios depends on the investment opportunities for the satellite portfolio. Today, we see a large imbalance between the supply of green assets and what is needed to finance the transition. Therefore, we can expect the allocation in the satellite portfolio to be small, but to increase in the future:

$$ \frac{\partial \alpha^{\text{satellite}}(t)}{\partial t} \geq 0 $$

The allocation $\alpha^{\text{satellite}}$ is then time-varying and must be revised periodically, e.g., once a year.

## 11.4 Exercise

### 11.4.1 Equity and bond portfolio optimization with green preferences

We consider an investment universe of 8 issuers. In the table below, we report the carbon emissions $\mathcal{C}\mathcal{E}_{i,j}$ (in ktCO$_2$e) of these companies and their revenues $Y_i$ (in $ bn)$, and we indicate in the last row whether the company belongs to sector $\mathcal{S}$ector$_1$ or $\mathcal{S}$ector$_2$:

<table><thead><tr><th>Issuer</th><th>#1</th><th>#2</th><th>#3</th><th>#4</th><th>#5</th><th>#6</th><th>#7</th><th>#8</th></tr></thead><tbody><tr><td><i>C</i><i>E</i><sub>i,1</sub></td><td>75</td><td>5 000</td><td>720</td><td>50</td><td>2 500</td><td>25</td><td>30 000</td><td>5</td></tr><tr><td><i>C</i><i>E</i><sub>i,2</sub></td><td>75</td><td>5 000</td><td>1 030</td><td>350</td><td>4 500</td><td>5</td><td>2000</td><td>64</td></tr><tr><td><i>C</i><i>E</i><sub>i,3</sub></td><td>24 000</td><td>15 000</td><td>1 210</td><td>550</td><td>500</td><td>187</td><td>30 000</td><td>199</td></tr><tr><td><i>Y</i><sub>i</sub></td><td>300</td><td>328</td><td>125</td><td>100</td><td>200</td><td>102</td><td>107</td><td>25</td></tr><tr><td><i>S</i>ector</td><td>1</td><td>2</td><td>1</td><td>1</td><td>2</td><td>1</td><td>2</td><td>2</td></tr></tbody></table>

The benchmark $b$ of this investment universe is defined as:

$$ b = (22\%, 19\%, 17\%, 13\%, 11\%, 8\%, 6\%, 4\%) $$

In what follows, we consider long-only portfolios.

1. We want to compute the carbon intensity of the benchmark.

(a) Compute the carbon intensities $\mathcal{C}\mathcal{I}_{i,j}$ of each company $i$ for the scopes 1, 2 and 3.

(b) Deduce the carbon intensities $\mathcal{C}\mathcal{I}_{i,j}$ of each company $i$ for the scopes 1 + 2 and 1 + 2 + 3.

(c) Deduce the weighted average carbon intensity (WACI) of the benchmark if we consider the scope 1 + 2 + 3.

(d) We assume that the market capitalization of the benchmark portfolio is equal to $10 bn and we invest $1 bn.
i. Deduce the market capitalization of each company (expressed in $ bn).

ii. Compute the ownership ratio for each asset (expressed in bps).

iii. Compute the carbon emissions of the benchmark portfolio<sup>43</sup> if we invest $1 bn and we consider the scope 1 + 2 + 3.

iv. Compare the (exact) carbon intensity of the benchmark portfolio with the WACI value obtained in Question 1.c.

2. We want to manage an equity portfolio with respect to the previous investment universe and reduce the weighted average carbon intensity of the benchmark by the rate $\mathcal{R}$. We assume that the volatility of the stocks is respectively equal to 22%, 20%, 25%, 18%, 40%, 23%, 13% and 29%. The correlation matrix between these stocks is given by:

$$ \rho = \begin{pmatrix} 100\% & & & & & & & \\ & 100\% & & & & & & \\ & & 75\% & 100\% & & & & \\ & & & 65\% & 80\% & 100\% & & \\ & & & & 50\% & 70\% & 85\% & 100\% \\ & & & & & 60\% & 70\% & 80\% \\ & & & & & & 50\% & 70\% & 75\% \\ & & & & & & & 60\% & 70\% & 75\% \\ & & & & & & & & 60\% & 100\% \end{pmatrix} $$

(a) Compute the covariance matrix $\Sigma$.

(b) Write the optimization problem if the objective function is to minimize the tracking error risk under the constraint of carbon intensity reduction.

(c) Give the QP formulation of the optimization problem.

(d) $\mathcal{R}$ is equal to 20%. Find the optimal portfolio if we target scope 1 + 2. What is the value of the tracking error volatility?

(e) Same question if $\mathcal{R}$ is equal to 30%, 50%, and 70%.

(f) We target scope 1 + 2 + 3. Find the optimal portfolio if $\mathcal{R}$ is equal to 20%, 30%, 50% and 70%. Give the value of the tracking error volatility for each optimized portfolio.

(g) Compare the optimal solutions obtained in Questions 2.e and 2.f.

3. We want to manage a bond portfolio with respect to the previous investment universe and reduce the weighted average carbon intensity of the benchmark by the rate $\mathcal{R}$. We use the scope 1 + 2 + 3. In the table below, we report the modified duration $MD_i$ and the duration-times-spread factor $DTS_i$ of each corporate bond $i$:

<table><thead><tr><th>Asset</th><th>#1</th><th>#2</th><th>#3</th><th>#4</th><th>#5</th><th>#6</th><th>#7</th><th>#8</th></tr></thead><tbody><tr><td>MD<sub>i</sub> (in years)</td><td>3.56</td><td>7.48</td><td>6.54</td><td>10.23</td><td>2.40</td><td>2.30</td><td>9.12</td><td>7.96</td></tr><tr><td>DTS<sub>i</sub> (in bps)</td><td>103</td><td>155</td><td>75</td><td>796</td><td>89</td><td>45</td><td>320</td><td>245</td></tr><tr><td><i>Sector</i></td><td>1</td><td>2</td><td>1</td><td>1</td><td>2</td><td>1</td><td>2</td><td>2</td></tr></tbody></table>

We remind that the active risk can be calculated using three functions. For the active share, we have:

$$ \mathcal{R}_{AS}(w | b) = \sigma_{AS}^2(w | b) = \sum_{i=1}^{n} (w_i - b_i)^2 $$

<footnote n="43"><sup>43</sup>We assume that the float percentage is equal to 100% for all the 8 companies.</footnote>
We also consider the MD-based tracking error risk:

$$ \mathcal{R}_{\text{MD}}(w | b) = \sigma_{\text{MD}}^2(w | b) = \sum_{j=1}^{n_{\text{Sector}}} \left( \sum_{i \in \text{Sector}_j} (w_i - b_i) \text{MD}_i \right)^2 $$

and the DTS-based tracking error risk:

$$ \mathcal{R}_{\text{DTS}}(w | b) = \sigma_{\text{DTS}}^2(w | b) = \sum_{j=1}^{n_{\text{Sector}}} \left( \sum_{i \in \text{Sector}_j} (w_i - b_i) \text{DTS}_i \right)^2 $$

Finally, we define the synthetic risk measure as a combination of AS, MD and DTS active risks:

$$ \mathcal{R}(w | b) = \varphi_{\text{AS}} \mathcal{R}_{\text{AS}}(w | b) + \varphi_{\text{MD}} \mathcal{R}_{\text{MD}}(w | b) + \varphi_{\text{DTS}} \mathcal{R}_{\text{DTS}}(w | b) $$

where $\varphi_{\text{AS}} \ge 0$, $\varphi_{\text{MD}} \ge 0$ and $\varphi_{\text{DTS}} \ge 0$ indicate the weight of each risk. In what follows, we use the following numerical values: $\varphi_{\text{AS}} = 100$, $\varphi_{\text{MD}} = 25$ and $\varphi_{\text{DTS}} = 1$. The reduction rate $\mathcal{R}$ of the weighted average carbon intensity is set to 50% for the scope 1 + 2 + 3.

(a) Compute the modified duration MD (b) and the duration-times-spread factor DTS (b) of the benchmark.

(b) Let $w_{\text{ew}}$ be the equally-weighted portfolio. Compute$^{44}$ MD ($w_{\text{ew}}$), DTS ($w_{\text{ew}}$), $\sigma_{\text{AS}}(w_{\text{ew}} | b)$, $\sigma_{\text{MD}}(w_{\text{ew}} | b)$ and $\sigma_{\text{DTS}}(w_{\text{ew}} | b)$.

(c) We consider the following optimization problem:

$$ \begin{aligned} w^* &= \arg \min \frac{1}{2} \mathcal{R}_{\text{AS}}(w | b) \\ \text{s.t.} \quad & \left\{ \begin{array}{l} \sum_{i=1}^{n} w_i = 1 \\ \text{MD}(w) = \text{MD}(b) \\ \text{DTS}(w) = \text{DTS}(b) \\ \mathcal{CI}(w) \le (1 - \mathcal{R}) \mathcal{CI}(b) \\ 0 \le w_i \le 1 \end{array} \right. \end{aligned} $$

Give the analytical value of the objective function. Find the optimal portfolio $w^*$. Compute MD ($w^*$), DTS ($w^*$), $\sigma_{\text{AS}}(w^* | b)$, $\sigma_{\text{MD}}(w^* | b)$ and $\sigma_{\text{DTS}}(w^* | b)$.

(d) We consider the following optimization problem:

$$ \begin{aligned} w^* &= \arg \min \frac{\varphi_{\text{AS}}}{2} \mathcal{R}_{\text{AS}}(w | b) + \frac{\varphi_{\text{MD}}}{2} \mathcal{R}_{\text{MD}}(w | b) \\ \text{s.t.} \quad & \left\{ \begin{array}{l} \sum_{i=1}^{n} w_i = 1 \\ \text{DTS}(w) = \text{DTS}(b) \\ \mathcal{CI}(w) \le (1 - \mathcal{R}) \mathcal{CI}(b) \\ 0 \le w_i \le 1 \end{array} \right. \end{aligned} $$

Give the analytical value of the objective function. Find the optimal portfolio $w^*$. Compute MD ($w^*$), DTS ($w^*$), $\sigma_{\text{AS}}(w^* | b)$, $\sigma_{\text{MD}}(w^* | b)$ and $\sigma_{\text{DTS}}(w^* | b)$.

<footnote>$^{44}$Precise the corresponding unit (years, bps or %) for each metric.</footnote>
(e) We consider the following optimization problem:

$$ w^* = \arg\min \frac{1}{2} \mathcal{R}(w | b) \\ \text{s.t.} \quad \left\{ \begin{array}{l} \sum_{i=1}^{n} w_i = 1 \\ \mathcal{CI}(w) \le (1 - \mathcal{R}) \mathcal{CI}(b) \\ 0 \le w_i \le 1 \end{array} \right. $$

Give the analytical value of the objective function. Find the optimal portfolio $w^*$. Compute MD($w^*$), DTS($w^*$), $\sigma_{AS}(w^* | b)$, $\sigma_{MD}(w^* | b)$ and $\sigma_{DTS}(w^* | b)$.

(f) Comment on the results obtained in Questions 3.c, 3.d and 3.e.

(g) How to find the previous solution of Question 3.e using a QP solver?

4. We consider a variant of Question 3 and assume that the synthetic risk measure is:

$$ \mathcal{D}(w | b) = \varphi_{\text{AS}}\mathcal{D}_{\text{AS}}(w | b) + \varphi_{\text{MD}}\mathcal{D}_{\text{MD}}(w | b) + \varphi_{\text{DTS}}\mathcal{D}_{\text{DTS}}(w | b) $$

where:

$$ \mathcal{D}_{\text{AS}}(w | b) = \frac{1}{2} \sum_{i=1}^{n} |w_i - b_i| $$

$$ \mathcal{D}_{\text{MD}}(w | b) = \sum_{j=1}^{n} \left| \sum_{i \in \mathcal{S}_{\text{sector}}_{j}} (w_i - b_i) \text{MD}_i \right| $$

$$ \mathcal{D}_{\text{DTS}}(w | b) = \sum_{j=1}^{n} \left| \sum_{i \in \mathcal{S}_{\text{sector}}_{j}} (w_i - b_i) \text{DTS}_i \right| $$

(a) Define the corresponding optimization problem when the objective is to minimize the active risk and reduce the carbon intensity of the benchmark by $\mathcal{R}$.

(b) Give the LP formulation of the optimization problem.

(c) Find the optimal portfolio when $\mathcal{R}$ is set to 50%. Compare the solution with this obtained in Question 3.e.

### 11.4.2 Monotonicity property of the order-statistic and naive approaches
Chapter 12

Physical Risk Modeling
## 12.1 General circulation model

12.1.1 Carbon cycle

12.1.2 Greenhouse gas chemistry

12.1.3 Horizontal and vertical heat transport

12.1.4 Oceans

12.1.5 Carbon sinks

## 12.2 Statistical modeling of climate hazards

12.2.1 Chronic risk

12.2.2 Acute risk
## 12.3 Geolocation

### 12.3.1 Climate hazard location

### 12.3.2 Asset location
12.4 Applications

12.4.1 Cyclones and hurricanes

12.4.2 Drought

12.4.3 Floods

12.4.4 Extreme heat

12.4.5 Water stress

12.4.6 Wildfire

12.5 Exercises
# Chapter 13

## Climate Stress Testing and Risk Management
# 13.1 Transmission channels

13.1.1 Direct and indirect transmission

13.1.2 Credit transmission channel

13.1.3 Market transmission channel

13.1.4 Systemic risk
13.2 Climate risk hedging

13.3 Climate value-at-risk

13.4 Climate stress testing

13.5 Exercises
<figure><img src="image_530.png" /></figure>
Conclusion
<figure><img src="image_531.png" /></figure>
# Appendices
<figure><img src="image_532.png" /></figure>
# Appendix A

# Technical Appendix

## A.1 Mathematical tools

### A.1.1 Linear algebra

#### Eigendecomposition

The value $\lambda$ is an eigenvalue of the $n \times n$ matrix $A$ if there exists a non-zero eigenvector $v$ such that we have $Av = \lambda v$. We denote $V$ the matrix composed of the $n$ eigenvectors. We have $AV = V\Lambda$ where $\Lambda = \text{diag}(\lambda_1, \dots, \lambda_n)$ is the diagonal matrix of eigenvalues. We finally obtain the eigendecomposition of the matrix $A$:

$$A = V\Lambda V^{-1}$$

If $A$ is an hermitian matrix<fnref n="1" />, then the matrix $V$ of eigenvectors is unitary. It follows that:

$$A = V\Lambda V^*$$

In particular, if $A$ is a symmetric real matrix, we obtain the following relationship<fnref n="2" />:

$$A = V\Lambda V^{\top}$$

#### Schur decomposition and matrix function

The Schur decomposition of the $n \times n$ matrix $A$ is equal to:

$$A = Q T Q^* \qquad (A.1)$$

where $Q$ is a unitary matrix and $T$ is an upper triangular matrix<fnref n="3" />. This decomposition is useful to calculate matrix functions.

Let us consider the matrix function in the space $\mathbb{M}$ of square matrices:

$$f: \mathbb{M} \to \mathbb{M} \\ A \mapsto B = f(A)$$

For instance, if $f(x) = \sqrt{x}$ and $A$ is positive, we can define the matrix $B$ such that:

$$BB^* = B^*B = A$$

<footnote n="1">¹The square matrix $A$ is hermitian if it is equal to its own conjugate transpose $A^*$, implying that we have $A_{i,j} = \text{conj} A_{j,i}$;</footnote>

<footnote n="2">²We have $A^\top = (V\Lambda V^{-1})^\top = (V^{-1})^\top \Lambda V^\top$. We deduce that $V^{-1} = V^\top$.</footnote>

<footnote n="3">³$Q$ and $T$ are also called the transformation matrix and the Schur form of $A$.</footnote>
B is called the square root of A and we note $B = A^{1/2}$. This matrix function generalizes the scalar-valued function to the set of matrices. Let us consider the following Taylor expansion:

$$f(x) = f(x_0) + (x - x_0) f'(x_0) + \frac{(x - x_0)^2}{2!} f''(x_0) + \dots$$

We can show that if the series converge for $|x - x_0| < \alpha$, then the matrix $f(A)$ defined by the following expression:

$$f(A) = f(x_0) + (A - x_0 I) f'(x_0) + \frac{(A - x_0 I)^2}{2!} f''(x_0) + \dots$$

converges to the matrix $B$ if $|A - x_0 I| < \alpha$ and we note $B = f(A)$. In the case of the exponential function, we have:

$$f(x) = e^x = \sum_{k=0}^{\infty} \frac{x^k}{k!}$$

We deduce that the exponential of the matrix $A$ is equal to:

$$B = e^A = \sum_{k=0}^{\infty} \frac{A^k}{k!}$$

In a similar way, the logarithm of $A$ is the matrix $B$ such that $e^B = A$ and we note $B = \ln A$.

Let $A$ and $B$ be two $n \times n$ square matrices. Using the Taylor expansion, Golub and Van Loan (2013) showed that $f(A^\top) = f(A)^\top$, $Af(A) = f(A)A$ and $f(B^{-1}AB) = B^{-1}f(A)B$. It follows that:

$$e^{A^\top} = (e^A)^\top$$

and:

$$e^{B^{-1}AB} = B^{-1}e^A B$$

If $AB = BA$, we can also prove that $Ae^B = e^B A$ and $e^{A+B} = e^A e^B = e^B e^A$.

**Remark 147** There are different ways to compute numerically $f(A)$. For transcendental functions, we have:

$$f(A) = Qf(T)Q^*$$

where $A = QTQ^*$ is the Schur decomposition of $A$. Because $T$ is an upper diagonal matrix, $f(T)$ is also a diagonal matrix whose elements can be calculated with Algorithm 9.1.1 of Golub and Van Loan (2013). This algorithm is reproduced below$^4$.

## Sherman-Morrison-Woodbury formula

Let $u$ and $v$ be two $n \times 1$ vectors and $A$ be an invertible $n \times n$ matrix. We can show that (Golub and Van Loan, 2013):

$$(A + uv^\top)^{-1} = A^{-1} - \frac{1}{1 + v^\top A^{-1} u} A^{-1} uv^\top A^{-1}$$

<footnote>$^4$For the exponential matrix, we may prefer to use the Pade approximation method, which is described in Algorithm 9.3.1 (scaling and squaring) of Golub and Van Loan (2013). See also the survey of Moler and Van Loan (2003).</footnote>
**Algorithm 2** Schur-Parlett matrix function f (A)

Compute the Schur decomposition $A = QTQ^*$

Initialize F to the matrix $0_{n \times n}$

for $i = 1:n$ do

$f_{i,i} \leftarrow f(t_{i,i})$

end for

for $p = 1:n-1$ do

for $i = 1:n-p$ do

$j \leftarrow i + p$

$s \leftarrow t_{i,j} (f_{j,j} - f_{i,i})$

for $k = i+1:j-1$ do

$s \leftarrow s + t_{i,k} f_{k,j} - f_{i,k} t_{k,j}$

end for

$f_{i,j} \leftarrow s / (t_{j,j} - t_{i,i})$

end for

end for

$B \leftarrow QFQ^*$

return B

Source: Golub and Van Loan (2013), page 519.

Batista and Karawia (2009) extended the SMW formula when the outer product is a sum:

$$ \left( A + \sum_{k=1}^{m} u_k v_k^{\top} \right)^{-1} = A^{-1} - A^{-1} U S^{-1} V^{\top} A^{-1} $$

where $U = (u_1 \cdots u_m)$ and $V = (v_1 \cdots v_m)$ are two $n \times m$ matrices, and $S = I_m + T$ and $T = (T_{i,j})$ are two $m \times m$ matrices where $T_{i,j} = v_i^{\top} A^{-1} u_j$. In the case $m = 2$, Roncalli et al. (2020) showed that the SMW formula becomes:

$$ (A + u_1 v_1^{\top} + u_2 v_2^{\top})^{-1} = A^{-1} - A^{-1} C A^{-1} $$

where:

$$ C = \frac{(I_n - u_2 v_2^{\top} A^{-1}) u_1 v_1^{\top} + (I_n - u_1 v_1^{\top} A^{-1}) u_2 v_2^{\top} + u_1 v_2^{\top} A^{-1} u_2 v_1^{\top} + u_2 v_1^{\top} A^{-1} u_1 v_2^{\top}}{1 + (v_1^{\top} A^{-1} u_1)(v_2^{\top} A^{-1} u_2) + v_1^{\top} A^{-1} u_1 + v_2^{\top} A^{-1} u_2 - (v_2^{\top} A^{-1} u_1)(v_1^{\top} A^{-1} u_2)} $$

### A.1.2 Optimization

#### Bisection algorithm

The simplest algorithm to find the roots of the equation $f(x) = 0$ is the bisection algorithm. We assume that the function $f$ is continuous on the interval $[a, b]$, and $f(a)$ and $f(b)$ have opposite signs: $f(a) f(b) < 0$. The underlying idea is to reduce the bracket by computing the midpoint $c = (a + b)/2$. If $f(a) f(c) < 0$, then the interval becomes $[a, c]$, otherwise it is equal to $[c, b]$. This step is repeated until the interval is sufficiently small. The absolute error is divided by 2 at each step so the method converges linearly. Let $c^*$ be a root such that $f(c^*) = 0$ and $c^* \in [a, b]$. Let $c_{(k)}$ be the midpoint at the $k^{th}$ iteration. Then the difference between $c_{(k)}$ and $c^*$ is bounded by $|c_{(k)} - c^*| \le 2^{-k} (b - a)$. For instance, after 10 iterations, we are sure that the interval has been divided by $2^{10} = 1024$. After 20 iterations, the interval is divided by $10^6$.
Linear programming

The canonical form of a linear program is:

$$
\begin{array}{ll}
x^* = \arg\max c^\top x & \\
\text{s.t.} \quad \left\{ \begin{array}{l} Sx \le T \\ x \ge \mathbf{0}_n \end{array} \right. &
\end{array}
$$

where $x$ is a $n \times 1$ vector, $c$ is a $n \times 1$ vector, $S$ is a $n_S \times n$ matrix, and $T$ is a $n_S \times 1$ vector.
A minimization linear programming (LP) problem can be cast into a maximization LP problem
because $\min c^\top x \equiv \max -c^\top x$. A maximization LP problem is in standard form if it is written as:

$$
\begin{array}{ll}
x^* = \arg\max c^\top x & \\
\text{s.t.} \quad \left\{ \begin{array}{l} Sx = T \\ x \ge \mathbf{0}_n \end{array} \right. &
\end{array}
$$

We notice that a canonical problem can always be put into a standard form by introducing the slack
variables $\xi$:

$$
Sx \le T \Leftrightarrow \{\xi \ge \mathbf{0}_{n_S} : Sx + I_{n_S}\xi = T\}
$$

We obtain the following LP problem:

$$
\begin{array}{ll}
y^* = \arg\max d^\top y & \\
\text{s.t.} \quad \left\{ \begin{array}{l} My = N \\ y \ge \mathbf{0}_n \end{array} \right. &
\end{array}
$$

where $d = (c, \mathbf{0}_{n_S})$, $y = (x, s)$, $M = (S \ I_{n_S})$ and $N = T$. The technique of slack variables<fnref n="5" /> can also be used when the constraint is $Sx \ge T$ and we have $Sx - \xi = T$. If we introduce lower and upper bounds ($x^- \le x \le x^+$), we can use the change of variable $y = x - x^-$. It follows that $y \ge \mathbf{0}_n$ and $I_n y \le x^+ - x^-$. Since we have $x = x^- + y$, we obtain $c^\top x = c^\top x^- + c^\top y$ and $Sx \le T \Leftrightarrow Sy \le T - Sx^-$. Finally, we deduce that:

$$
\left\{
\begin{array}{l}
x^* = \arg\max c^\top x \\
\text{s.t.} \quad \left\{ \begin{array}{l} Sx \le T \\ x \ge \mathbf{0}_n \end{array} \right.
\end{array}
\right.
\iff
\left\{
\begin{array}{l}
y^* = \arg\max c^\top y \\
\text{s.t.} \quad \left\{ \begin{array}{l} \left( \begin{array}{c} S \\ I_n \end{array} \right) y \le \left( \begin{array}{c} T - Sx^- \\ x^+ - x^- \end{array} \right) \\ y \ge \mathbf{0}_n \end{array} \right.
\end{array}
\right.
$$

If some variables $x_i \in \mathbb{R}$, we can write $x_i = y_i - z_i$ where $y_i \ge 0$ and $z_i \ge 0$.

The previous analysis shows that we can always transform a general LP problem into a canonical
or standard form. This is why most numerical packages consider the following general formulation:

$$
\begin{equation}
\begin{split}
x^* = \arg\min c^\top x & \qquad \text{(A.2)} \\
\text{s.t.} \quad \left\{
\begin{array}{l}
Ax = B \\
Cx \leq D \\
x^- \leq x \leq x^+
\end{array}
\right.
&
\end{split}
\end{equation}
$$

where $x$ is a $n \times 1$ vector, $c$ is a $n \times 1$ vector, $A$ is a $n_A \times n$ matrix, $B$ is a $n_A \times 1$ vector, $C$ is a $n_C \times n$ matrix, $D$ is a $n_C \times 1$ vector, and $x^-$ and $x^+$ are two $n \times 1$ vectors. If $n_A = 0$, there is no equality constraints. Similarly, there is no inequality constraints if $n_C = 0$. If there is no lower

<footnote n="5">⁵When we rewrite an inequality by subtracting new variables, these last ones are also called surplus variables.</footnote>
bounds or/and upper bounds, this implies that $x^{-} = -\infty \cdot \mathbf{1}_{n}$ and $x^{+} = \infty \cdot \mathbf{1}_{n}$. A solution $x$ is said to be feasible if it satisfies all the constraints of the linear program (A.2). The set of feasible solutions is called the feasible space $\Omega$:

$$ \Omega = \{x \in \mathbb{R}^n : Ax = B, Cx \leq D, x^{-} \leq x \leq x^{+}\} $$

Since $Ax = B$ and $Cx \leq D$ define two convex sets (a hyperplane and a closed half-space), $\Omega$ is a convex set. If $\Omega$ is not empty and $c^{\top}x$ is bounded below, then there is an optimal solution.

**Example 65** We consider the following optimization problem:

$$ \begin{array}{ll} (x_1^*, x_2^*, x_3^*) = & \arg\min x_1 + x_2 + x_3 \\ \text{s.t.} \quad & \begin{cases} x_1 + x_2 + x_3 = 1 \\ 3x_1 + x_2 + 5x_3 \geq 2 \\ 2x_1 + x_2 + 3x_3 \leq s \\ 0 \leq x_i \leq 1 \end{cases} \end{array} $$

If $s \geq 1.5$, the optimal solution is $(0, 0.75, 0.25)$. If $s < 1.5$, $\Omega = \emptyset$ and there is no solution.

Let us consider the standard program: $x^* = \arg\max f(x) = c^\top x$ subject to $Ax = b$. The Lagrange function is:

$$ \mathcal{L}(x; \lambda) = c^{\top}x - \lambda^{\top}(Ax - b) $$

The first-order condition is:

$$ \left\{ \begin{array}{l} \partial_x \mathcal{L}(x; \lambda) = c - A^{\top}\lambda = \mathbf{0}_n \\ \partial_{\lambda} \mathcal{L}(x; \lambda) = Ax - b = \mathbf{0}_{n_A} \end{array} \right. $$

The optimal solution $(x^*, \lambda^*)$ satisfies the following system of equations: $Ax^* = b$ and $A^\top \lambda^* = c$. We also have $f(x^*) = \mathcal{L}(x^*, \lambda^*) = c^\top x^* - \lambda^* \top (Ax^* - b)$ and:

$$ \mathcal{L}(x^*, \lambda^*) = \sum_{i=1}^{n} c_i x_i^* - \sum_{j=1}^{n_A} \lambda_j^* \sum_{i=1}^{n} A_{j,i} x_i^* + \sum_{j=1}^{n_A} \lambda_j^* b_j $$

Let us assume that we increase $b$ by a vector $\varepsilon \geq \mathbf{0}_{n_A}$. The Lagrange function becomes:

$$ \begin{align*} \mathcal{L}_{\varepsilon}(x^*, \lambda^*) &= \sum_{i=1}^{n} c_i x_i^* - \sum_{j=1}^{n_A} \lambda_j^* \sum_{i=1}^{n} A_{j,i} x_i^* + \sum_{j=1}^{n_A} \lambda_j^* (b_j + \varepsilon_j) \\ &= f(x^*) + \sum_{j=1}^{n_A} \lambda_j^* \varepsilon_j \end{align*} $$

If we relax the $j^{th}$ constraint by a small value of $\varepsilon_j$, we change the objective function by $\lambda_j^* \varepsilon_j$. If $\lambda_j^* = 0$, then modifying locally the constraint has no impact. If $\lambda_j^* < 0$, we have $f_{\varepsilon}(x^*) < f(x^*)$, implying that we have no incentive to change the value $b_j$ of the constraint, otherwise we will reduce the objective function. If $\lambda_j^* > 0$, we have $f_{\varepsilon}(x^*) > f(x^*)$, meaning that we can improve the objective function. From an economic point of view, a linear program is usually a profit maximization, a production maximization or a cost minimization. In this case, we can ask ourselves what price $p_j$ we are willing to pay to increase the right-hand side value $b_j$ by one unit? By taking $\varepsilon_j = 1$, we obtain $p_j \leq \lambda_j^*$ and $\lambda_j^*$ is called the shadow or dual price. This explains why duality plays an important role in economics.
**Example 66** We consider the following optimization problem:

$$ 
\begin{array}{l} 
x^* = \arg \max 2x_2 + 3x_3 - x_1 \\ 
\text{s.t.} \quad \begin{cases} x_1 + x_2 + x_3 = 1 \\ 2x_1 + 5x_2 + 6x_3 = 0.2 \end{cases} 
\end{array} 
$$

The optimal solution is $x^* = (1.45, 0, -0.45)$ and $\lambda^* = (3, -1)$, and the objective function is equal to $f(x^*) = 2.80$. The shadow prices are then equal to +3 for the first constraint and -1 for the second constraint. We do not want to increase the second constraint by $\varepsilon_2$ because it reduces the objective function by $\varepsilon_2$. Indeed, if $\varepsilon_2 = 0.04$, the new solution is $x^* = (1.44, 0, -0.44)$ and we have $f(x^*) = 2.76$. On the contrary, if we increase the first constraint by $\varepsilon_1$, we improve the objective function by $3\varepsilon_1$. For instance, if $\varepsilon_1 = 0.1$, the new solution is $x^* = (1.6, 0, -0.5)$ and we have $f(x^*) = 3.1$. Let us combine the two effects. If $(\varepsilon_1, \varepsilon_2) = (0.1, 0.04)$, we have $x^* = (1.59, 0, -0.49)$ and $f(x^*) = 3.06$. We verify that:

$$ 
\begin{aligned} 
f(x^*) + \lambda_1^* \varepsilon_1 + \lambda_2^* \varepsilon_2 &= 2.80 + 3 \times 0.1 - 1 \times 0.04 \\ 
&= 3.06 \\ 
&= f(x^*) 
\end{aligned} 
$$

If the primal problem is maximizing $c^\top x$ subject to $Ax \le b$ and $x \ge 0_n$, then the dual problem is defined as:

$$ 
\begin{array}{l} 
x^* = \arg \min b^\top y \\ 
\text{s.t.} \quad \begin{cases} A^\top y \ge c \\ y \ge 0_m \end{cases} 
\end{array} 
\qquad (A.3) 
$$

We notice that:

* The dual problem has a linear programming form;

* The maximization problem and the $\le$ inequality type are changed into a minimization problem and a $\ge$ inequality type;

* Each dual variable $y_j$ may be assigned to a corresponding primal constraint $\sum_{i=1}^n A_{j,i}x_i \le b_j$;

* The roles of the vectors $c$ and $b$ are switched: $c$ becomes the right-hand side of the constraints and $b$ the coefficients of the objective function;

* The matrix $A$ of the inequality constraints is transposed.

Let $\Omega_x = \{x : Ax \le b, x \ge 0_n\}$ and $\Omega_y = \{y : A^\top y \ge c, y \ge 0_m\}$ be the primal and dual feasible spaces. The weak duality theorem states that$^6$:

$$ \forall x \in \Omega_x, y \in \Omega_y : c^\top x \le b^\top y $$

This means that the objective value of the dual problem is an upper bound on the objective value of the primal problem, and vice versa. Using this property, we can also formulate the strong duality theorem: A primal feasible solution $x^*$ is optimal if and only if there exists a dual feasible solution $y^*$ such that $c^\top x^* = b^\top y^*$. We also deduce that $y^*$ is the optimal solution of the dual LP problem.

<footnote>$^6$If $x \ge 0_n$ and $\alpha \ge \beta$, we have $\alpha^\top x \ge \beta^\top x$ because $\alpha_i x_i \ge \beta_i x_i$ and $\sum_{i=1}^n \alpha_i x_i \ge \sum_{i=1}^n \beta_i x_i$. Using the constraints $Ax \le b$ and $A^\top y \ge c$, we deduce that $c^\top x = x^\top c \le x^\top A^\top y = (Ax)^\top y \le b^\top y$.</footnote>
**Example 67** We consider a classical ressource allocation problem. We would like to produce 3 goods, whose market prices are $3, $5 and $7. The production of these goods requires some ressources, such as capital, labor, raw materials and land. The available ressources are respectively equal to 100, 50, 80 and 200. For the first good, the production of one unit needs 5 units of capital, 3 units of labor, 1 unit of raw materials and 1 unit of land. For the second good, the intermediary consumptions are 10, 2, 10 and 2, while they are equal to 4, 7, 8 and 4 for the third good.

The problem is to maximize the revenue under allocation constraints. By definition, the revenue is the product of sold quantities and market prices:

$$Revenue = 3x_2 + 5x_3 + 7x_1$$

If we consider the first ressource constraint, we have 100 units of capital that we can allocate between the three goods in the following way. We need 5 units of capital to produce one unit of the first good, 10 units of capital to produce one unit of the second good and 4 units of capital to produce one unit of the third good. The optimization problem is then defined as follows:

$$x^* = \arg\max 3x_2 + 5x_3 + 7x_1 \\ \text{s.t.} \quad \begin{cases} 5x_1 + 10x_2 + 4x_3 \le 100 \\ 3x_1 + 2x_2 + 7x_3 \le 50 \\ x_1 + 10x_2 + 8x_3 \le 80 \\ x_1 + 2x_2 + 4x_3 \le 200 \\ x_1 \ge 0, x_2 \ge 0, x_3 \ge 0 \end{cases}$$

The optimal solution is $x_1^* = 7.4390$, $x_2^* = 5.3049$ and $x_3^* = 2.4390$. The optimal value of the revenue is $65.91. The computation of the Lagrange multipliers gives $\lambda_1^* = 0.1555$, $\lambda_2^* = 0.6707$, $\lambda_3^* = 0.2104$ and $\lambda_4^* = 0$. While the primal LP was a resource allocation problem, the dual LP is a resource valuation problem from an economic viewpoint. Indeed, $y_i$ is the implicit cost per unit associated with each resource $i$ and the objective is then to minimize the total cost of the resources:

$$Cost = 100y_1 + 50y_2 + 80y_3 + 200y_4$$

Moreover, if we consider the first good, its market price $p_1$ is $3, and its internal value to produce it is equal to $v_1 = 5y_1 + 3y_2 + y_3 + y_4$. Therefore, we have the following constraint $v_1 \ge p_1$. Finally, we deduce that the resource valuation problem is to minimize the total cost under the constraints $v_1 \ge p_1$, $v_2 \ge p_2$ and $v_3 \ge p_3$. Therefore, the dual problem is:

$$y^* = \arg\min 100y_1 + 50y_2 + 80y_3 + 200y_4 \\ \text{s.t.} \quad \begin{cases} 5y_1 + 3y_2 + y_3 + y_4 \ge 3 \\ 10y_1 + 2y_2 + 10y_3 + 2y_4 \ge 5 \\ 4y_1 + 7y_2 + 8y_3 + 4y_4 \ge 7 \\ y_1 \ge 0, y_2 \ge 0, y_3 \ge 0, y_4 \ge 0 \end{cases}$$

The optimal solution is $y_1^* = 0.1555$, $y_2^* = 0.6707$, $y_3^* = 0.2104$ and $y_4^* = 0$. The optimal value of the dual objective function is $65.91$. The computation of the Lagrange multipliers gives $\vartheta_1^* = 7.4390$, $\vartheta_2^* = 5.3049$ and $\vartheta_3^* = 2.4390$. We can make the following observations. First, the strong duality property is satisfied because $Revenue(x^*) = Cost(y^*) = 65.91$. Second, we verify that $y^* = \lambda^*$ and $x^* = \vartheta^*$. The optimal solution of one problem is equal to the Lagrange coefficients of the second problem. Therefore, we interpret $y^*$ as the shadow prices of the resources. Third, a shadow price $y_j^*$ is equal to zero when the $j^{th}$ constraint is non-binding. Otherwise, it is strictly positive. If we
consider the previous example, the shadow price of the fourth resource is equal to zero because this resource is not scarce. In particular, we do not need to buy supplementary land to increase the revenue. Therefore, we are ready to pay zero for additional land. On the contrary, we can pay the supplementary unit up to $0.1555 to buy more capital. If the price of capital is greater than $0.1555, the increase of revenue will be lower than the increase of cost and we have no interest to increase the production capacity.

## Quadratic programming

A quadratic programming (QP) problem is an optimization problem with a quadratic objective function and linear inequality constraints:

$$ 
\begin{aligned}
x^* &= \arg \min \frac{1}{2} x^\top Q x - x^\top R \\
\text{s.t.} \quad & S x \le T
\end{aligned}
\qquad (A.4)
$$

where $x$ is a $n \times 1$ vector, $Q$ is a $n \times n$ matrix and $R$ is a $n \times 1$ vector. We note that the system of constraints $Sx \le T$ allows specifying linear equality constraints$^7$ $Ax = B$ or weight constraints $x^{-} \le x \le x^{+}$. Most numerical packages then consider the following formulation:

$$ 
\begin{aligned}
x^* &= \arg \min \frac{1}{2} x^\top Q x - x^\top R \\
\text{s.t.} \quad & \begin{cases} Ax = B \\ Cx \le D \\ x^{-} \le x \le x^{+} \end{cases}
\end{aligned}
\qquad (A.5)
$$

because the problem (A.5) is equivalent to the canonical problem (A.4) with the following system of linear inequalities:

$$ 
\begin{bmatrix} -A \\ A \\ C \\ -I_n \\ I_n \end{bmatrix} \quad x \le \begin{bmatrix} -B \\ B \\ D \\ -x^{-} \\ x^{+} \end{bmatrix}
$$

If the space $\Omega$ defined by $Sx \le T$ is non-empty and if $Q$ is a symmetric positive definite matrix, the solution exists because the function $f(x) = \frac{1}{2}x^\top Q x - x^\top R$ is convex. In the general case where $Q$ is a square matrix, the solution may not exist.

The Lagrange function is also:

$$ \mathcal{L}(x; \lambda) = \frac{1}{2} x^\top Q x - x^\top R + \lambda^\top (Sx - T) $$

We deduce that the dual problem is defined by:

$$ 
\begin{aligned}
\lambda^* &= \arg \max \left\{ \inf_x \mathcal{L}(x; \lambda) \right\} \\
\text{s.t.} \quad & \lambda \ge 0
\end{aligned}
$$

<footnote>$^7$This is equivalent to imposing that $Ax \ge B$ and $Ax \le B$.</footnote>
We note that $\partial_x \mathcal{L}(x; \lambda) = Qx - R + S^\top \lambda$. The solution to the problem $\partial_x \mathcal{L}(x; \lambda) = 0$ is then $x = Q^{-1}(R - S^\top \lambda)$. We obtain:

$$
\begin{align*}
\inf_{x} \mathcal{L}(x; \lambda) &= \frac{1}{2} (R^{\top} - \lambda^{\top} S) Q^{-1} (R - S^{\top} \lambda) - (R^{\top} - \lambda^{\top} S) Q^{-1} R + \\
&\quad \lambda^{\top} (SQ^{-1} (R - S^{\top} \lambda) - T) \\
&= \frac{1}{2} R^{\top} Q^{-1} R - \lambda^{\top} SQ^{-1} R + \frac{1}{2} \lambda^{\top} SQ^{-1} S^{\top} \lambda - R^{\top} Q^{-1} R + \\
&\quad 2 \lambda^{\top} SQ^{-1} R - \lambda^{\top} SQ^{-1} S^{\top} \lambda - \lambda^{\top} T \\
&= -\frac{1}{2} \lambda^{\top} SQ^{-1} S^{\top} \lambda + \lambda^{\top} (SQ^{-1} R - T) - \frac{1}{2} R^{\top} Q^{-1} R
\end{align*}
$$

The dual program is another quadratic program:

$$
\begin{equation}
\begin{aligned}
\lambda^* &= \arg \min \frac{1}{2} \lambda^\top \bar{Q} \lambda - \lambda^\top \bar{R} && (A.6) \\
\text{s.t.} \quad \lambda \ge 0
\end{aligned}
\end{equation}
$$

where $\bar{Q} = SQ^{-1}S^\top$ and $\bar{R} = SQ^{-1}R - T$.

Alternating direction method of multipliers

The alternating direction method of multipliers (ADMM) is an algorithm introduced by Gabay and
Mercier (1976) to solve optimization problems that can be expressed as:

$$
\begin{equation}
\begin{aligned}
\{x^*, y^*\} &= \arg \min_{(x,y)} f_x(x) + f_y(y) && (A.7) \\
\text{s.t.} \quad &Ax + By = c
\end{aligned}
\end{equation}
$$

where $A \in \mathbb{R}^{p \times n}$, $B \in \mathbb{R}^{p \times m}$, $c \in \mathbb{R}^p$, and the functions $f_x : \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$ and $f_y : \mathbb{R}^m \to \mathbb{R} \cup \{+\infty\}$ are proper closed convex functions. Boyd et al. (2011) showed that the ADMM algorithm consists of the following three steps:

1. The x-update is:

$$
x^{(k+1)} = \arg \min_x \left\{ f_x^{(k+1)}(x) := f_x(x) + \frac{\varphi}{2} \left\| Ax + By^{(k)} - c + u^{(k)} \right\|_2^2 \right\} \quad (A.8)
$$

2. The y-update is:

$$
y^{(k+1)} = \arg \min_y \left\{ f_y^{(k+1)}(y) := f_y(y) + \frac{\varphi}{2} \left\| Ax^{(k+1)} + By - c + u^{(k)} \right\|_2^2 \right\} \quad (A.9)
$$

3. The u-update is:

$$
u^{(k+1)} = u^{(k)} + (Ax^{(k+1)} + By^{(k+1)} - c) \quad (A.10)
$$

In this approach, $u^{(k)}$ is the dual variable of the primal residual $r = Ax + By - c$ and $\varphi$ is the $\mathcal{L}_2$-norm penalty variable. The parameter $\varphi$ may be constant or it may change at each iteration. The ADMM algorithm benefits from the dual ascent principle and the multipliers method. The difference with the latter is that the $x$- and $y$-updates are performed alternately. This makes it more flexible, since the updates are equivalent to computing proximal operators for $f_x$ and $f_y$ independently. In practice, ADMM may be slow to converge at high accuracy, but is fast to converge at modest accuracy. Therefore, ADMM is a good candidate for solving large-scale machine learning problems where high accuracy does not necessarily lead to a better solution.
### A.1.3 Quadratic form

#### Definition

In mathematics, a quadratic form is a polynomial with terms all of degree two:

$$QF(x_1, \dots, x_n) = \sum_{i=1}^{n} \sum_{j=1}^{n} a_{i,j} x_i x_j = x^\top A x$$

where $A = (a_{i,j})$ is a $n \times n$ matrix. Since $x^\top A x$ is a scalar, we have $(x^\top A x)^\top = x^\top A^\top x$ and:

$$ 
\begin{aligned}
QF(x_1, \dots, x_n) &= \frac{1}{2} (x^\top A x + x^\top A^\top x) \\
&= \frac{1}{2} x^\top (A + A^\top) x \\
&= \frac{1}{2} x^\top Q x
\end{aligned}
$$

where $Q = A + A^\top$ is a symmetric matrix. We use the notation $QF(x; Q) = \frac{1}{2} x^\top Q x$ to define a canonical quadratic form. A generalized quadratic form is a polynomial including terms with degrees one and zero<fnref n="8" />:

$$QF(x; Q, R, c) = \frac{1}{2} x^\top Q x - x^\top R + c$$

#### Main properties

We list here some properties that are helpful when considering portfolio optimization:

* The multiplication of a quadratic form by a scalar $\varphi \ge 0$ remains a quadratic form:

$$\varphi \cdot QF(w; Q, R, c) = QF(w; \varphi Q, \varphi R, \varphi c)$$

* The sum of two quadratic forms is a quadratic form:

$$QF(x; Q_1, R_1, c_1) + QF(x; Q_2, R_2, c_2) = QF(x; Q_1 + Q_2, R_1 + R_2, c_1 + c_2)$$

* A quadratic form applied to the difference vector $x - y$ is a quadratic form in $x$:

$$ 
\begin{aligned}
QF(x - y; Q, R, c) &= \frac{1}{2} (x - y)^\top Q (x - y) - (x - y)^\top R + c \\
&= \frac{1}{2} (x^\top Q x - 2x^\top Qy + y^\top Qy) - x^\top R + y^\top R + c \\
&= \frac{1}{2} x^\top Q x - x^\top (R + Qy) + \frac{1}{2} y^\top Qy + y^\top R + c \\
&= QF \left( x; Q, R + Qy, \frac{1}{2} y^\top Qy + y^\top R + c \right)
\end{aligned}
$$

and a quadratic form in $y$:

$$ 
\begin{aligned}
QF(x - y; Q, R, c) &= \frac{1}{2} y^\top Q y - y^\top (Qx - R) + \frac{1}{2} x^\top Q x - x^\top R + c \\
&= QF \left( y; Q, Qx - R, \frac{1}{2} x^\top Q x - x^\top R + c \right)
\end{aligned}
$$

<footnote n="8">⁸We add the one-degree polynomial $-x^\top R$ instead of $+x^\top R$ in order to match the canonical form of a QP problem.</footnote>
* When $Q$ is a diagonal matrix, the quadratic form reduces to:

$$Q\mathcal{F}(x; Q, R, c) = \frac{1}{2} \sum_{i=1}^{n} q_{i,i} x_{i}^{2} - \sum_{i=1}^{n} r_{i} x_{i} + c$$

We deduce that:

$$\frac{1}{2} \sum_{i=1}^{n} q_{i} x_{i}^{2} = Q\mathcal{F}(x; \mathcal{D}(q), \mathbf{0}_{n}, 0)$$

where $q = (q_1, \dots, q_n)$ is a $n \times 1$ vector and $\mathcal{D}(q) = \text{diag}(q)$.

* Using the previous properties, we have:

$$\frac{1}{2} \sum_{i=1}^{n} q_i (x_i - y_i)^2 = Q\mathcal{F}(x - y; \mathcal{D}(q), \mathbf{0}_n, 0) \\ = Q\mathcal{F}\left(x; \mathcal{D}(q), \mathcal{D}(q) y, \frac{1}{2} y^\top \mathcal{D}(q) y\right)$$

* The square of the weighted sum is also a quadratic form:

$$\frac{1}{2} \left( \sum_{i=1}^{n} q_i x_i \right)^2 = \frac{1}{2} \left( x^\top q \right)^2 = \frac{1}{2} \left( x^\top q \right) \left( x^\top q \right)^\top = \frac{1}{2} x^\top q q^\top x^\top = Q\mathcal{F}(x; T(q), \mathbf{0}_n, 0)$$

where $T(q) = qq^\top$.

* We deduce that:

$$\frac{1}{2} \left( \sum_{i=1}^{n} q_i (x_i - y_i) \right)^2 = Q\mathcal{F}(x - y; T(q), \mathbf{0}_n, 0) \\ = Q\mathcal{F}\left(x; T(q), T(q) y, \frac{1}{2} y^\top T(q) y\right)$$

* The previous results can be extended when we consider partial sums:

$$\sum_{i \in \Omega} q_i x_i = \sum_{i=1}^{n} 1\{i \in \Omega\} \cdot q_i x_i = \sum_{i=1}^{n} (\omega_i q_i) x_i = \sum_{i=1}^{n} \tilde{q}_i x_i$$

where $\omega_i = 1\{i \in \Omega\}$ and $\tilde{q}_i = \omega_i q_i$. In a matrix form, we have:

$$\sum_{i \in \Omega} q_i x_i = (\omega \circ q)^\top x$$

where $\omega = (\omega_1, \dots, \omega_n)$. We deduce that:

$$\left\{ \begin{array}{l} \frac{1}{2} \sum_{i \in \Omega} q_i x_i^2 = Q\mathcal{F}(x; \mathcal{D}(\omega \circ q), \mathbf{0}_n, 0) \\ \\ \frac{1}{2} \sum_{i \in \Omega} q_i (x_i - y_i)^2 = Q\mathcal{F}\left(x; \mathcal{D}(\omega \circ q), \mathcal{D}(\omega \circ q) y, \frac{1}{2} y^\top \mathcal{D}(\omega \circ q) y\right) \\ \\ \frac{1}{2} \left(\sum_{i \in \Omega} q_i x_i\right)^2 = Q\mathcal{F}(x; T(\omega \circ q), \mathbf{0}_n, 0) \\ \\ \frac{1}{2} \left(\sum_{i \in \Omega} q_i (x_i - y_i)\right)^2 = Q\mathcal{F}\left(x; T(\omega \circ q), T(\omega \circ q) y, \frac{1}{2} y^\top T(\omega \circ q) y\right) \end{array} \right.$$

We notice that $\mathcal{D}(\omega \circ q) = \text{diag}(\omega \circ q) = \mathcal{D}(\omega)\mathcal{D}(q)$ and $T(\omega \circ q) = (\omega \circ q)(\omega \circ q)^\top = (\omega\omega^\top) \circ qq^\top = T(\omega) \circ T(q)$.
### A.1.4 Nonnegative matrix

#### Definition

$A = (a_{i,j})$ is a nonnegative matrix if all the elements are greater than or equal to zero: $a_{i,j} \ge 0$. Here are some special cases:

* A positive matrix is one where all elements are strictly greater than zero.

* A stochastic matrix is a square nonnegative matrix whose rows or columns are probability vectors<fnref n="9" />.

- A right/row stochastic matrix is a stochastic matrix where each row sums to 1: $\sum_{j=1}^{n} a_{i,j} = 1$.

- A left/column stochastic matrix is a stochastic matrix where each column sums to 1: $\sum_{i=1}^{n} a_{i,j} = 1$.

- A doubly stochastic matrix is a right and left stochastic matrix.

* A substochastic matrix is a square nonnegative matrix whose rows or columns add up to at most 1.

- A is a right/row substochastic matrix if $\sum_{j=1}^{n} a_{i,j} \le 1$.

- A is a left/column substochastic matrix if $\sum_{i=1}^{n} a_{i,j} \le 1$.

- A is a doubly substochastic matrix if it is a right and left substochastic matrix.

#### Ordering properties

$A$ is greater than or equal to $B$ if and only if $A_{i,j} \ge B_{i,j}$ for all $i, j$. We use the notation $A \succeq B$. By definition, a square nonnegative matrix $A$ satisfies $A \succeq 0_{n,n}$. If $A$ is a stochastic matrix, we have $0_{n,n} \preceq A \preceq 1_{n,n}$. Let $A, B, C$ and $D$ be square nonnegative matrices. We can show that:

$$ (NN1) \quad A \succeq B \Rightarrow AC \succeq BC; $$

$$ (NN2) \quad A \succeq B \wedge C \succeq D \Rightarrow AC \succeq BD; $$

$$ (NN3) \quad A \succeq B \wedge C \succeq D \Rightarrow A + C \succeq B + D; $$

$$ (NN4) \quad A \succeq B \wedge k \ge 1 \Rightarrow A^k \succeq B^k; $$

The proofs can be found in Desnos et al. (2023).

<footnote n="9">⁹This means that each entry is a real number between 0 and 1 and the sum of entries is equal to 1.</footnote>
## A.2 Statistical and probability analysis

### A.2.1 Probability distributions

#### The Bernoulli distribution

The Bernoulli random variable $X$ takes the value 1 with success probability of $p$ and the value 0 with failure probability of $q = 1 - p$. We note $X \sim \mathcal{B}(p)$. The probability mass function may also be expressed as follows:

$$ \Pr\{X = k\} = p^k (1-p)^{1-k} \quad \text{with } k = 0, 1 $$

We have $\mathbb{E}[X] = p$ and $\text{var}(X) = p(1-p)$.

#### The binomial distribution

The binomial random variable $X$ is the sum of $n$ independent Bernoulli random variables with the same probability of success $p$:

$$ X = \sum_{i=1}^{n} \mathcal{B}_i(p) $$

We note $X \sim \mathcal{B}(n, p)$. The probability mass function is equal to:

$$ \Pr\{X = k\} = \binom{n}{k} p^k (1-p)^{n-k} \quad \text{with } k = 0, 1, \dots, n $$

We have $\mathbb{E}[X] = np$ and $\text{var}(X) = np(1-p)$.

#### The geometric distribution

The geometric random variable $X$ is the number of Bernoulli trials needed to get one success. We note $X \sim \mathcal{G}(p)$. The probability mass function is equal to:

$$ \Pr\{X = k\} = (1-p)^{k-1} p \quad \text{with } k \in \mathbb{N}^* $$

We have $\mathbb{E}[X] = 1/p$ and $\text{var}(X) = (1-p)/p^2$.

**Remark 148** If we define $X$ as the number of failures before the first success, we have $\Pr\{X = k\} = (1-p)^k p$ with $k \in \mathbb{N}$, $\mathbb{E}[X] = (1-p)/p$ and $\text{var}(X) = (1-p)/p^2$.

#### The Poisson distribution

The Poisson random variable $X$ is the number of times an event occurs in the unit interval of time. We note $X \sim \mathcal{P}(\lambda)$ where $\lambda$ is the parameter of the Poisson distribution. The probability mass function is equal to:

$$ \Pr\{X = k\} = \frac{\lambda^k e^{-\lambda}}{k!} \quad \text{with } k \in \mathbb{N} $$

We have $\mathbb{E}[X] = \text{var}(X) = \lambda$. The parameter $\lambda$ is then the expected number of events occurring in the unit interval of time.
## The 1- and n-diversity distributions

Let $\pi_n \in \mathbb{R}_+^n$ such that $1^\top \pi_n = 1$. $\pi_n$ is then a probability distribution. The probability distribution $\pi_n$ is perfectly concentrated if there exists one observation $i_0$ such that $\pi_{n,i_0} = 1$ and $\pi_{n,i} = 0$ if $i \neq i_0$. When $n$ tends to $+\infty$, the limit distribution is noted $\pi_\infty$. On the opposite, the probability distribution $\pi_n^+$ such that $\pi_{n,i}^+ = 1/n$ for all $i = 1, \dots, n$ has no concentration. $\pi_n^-$ and $\pi_n^+$ are respectively called 1- and n-diversity distributions.

## The negative binomial distribution

### Probability distribution!Negative binomial distribution

The negative binomial distribution is a probability distribution used for modeling the frequency of an event. We denote it as $X \sim \mathcal{NB}(r, p)$, where $r > 0$ and $p \in [0, 1]$. The probability mass function is equal to:

$$ \Pr\{X = k\} = \binom{r+k-1}{k} (1-p)^r p^k \quad \text{with } k \in \mathbb{N} $$

We have $\mathbb{E}[X] = pr/(1-p)$ and $\text{var}(X) = pr/(1-p)^2$.

## The Poisson binomial distribution

The Poisson binomial distribution is the discrete probability distribution of a sum of independent Bernoulli trials that are not identically distributed $^{10}$:

$$ X = \sum_{i=1}^{n} \mathcal{B}(p_i) \sim \mathcal{PB}(p_1, \dots, p_n) $$

The probability mass function is given by:

$$ \Pr\{X = k\} = \sum_{\mathcal{E} \in \mathcal{F}_k} \prod_{i \in \mathcal{E}} p_i \prod_{j \in \mathcal{E}^c} (1-p_j) $$

where $\mathcal{F}_k$ is the set of all subsets of $k$ integers that can be selected from $\{1, \dots, n\}$ and $\mathcal{E}^c$ denotes the complement of the subset $\mathcal{E}$. The probability mass function can also be expressed as the convolution of the vectors $\{u_0, u_1, \dots, u_n\}$ where $u_0 = 1$ and $u_i = (p_i, 1-p_i)$. To calculate $\Pr\{X = k\}$, we deduce that the vector of discrete probabilities is the inverse Fourier transform of the vector $v_n$ where:

$$ v_j = \left( 1 + p_j + p_j \exp \left( -\frac{2\pi i}{n+1} \right) \right) v_{j-1} $$

and $v_0 = \mathbf{1}_{n+1}$. We have $\mathbb{E}[X] = \sum_{i=1}^n p_i$ and $\text{var}(X) = \sum_{i=1}^n p_i (1-p_i)$.

## The uniform distribution

$X$ is a uniform random variable $\mathcal{U}_{[a,b]}$ if the density function is $f(x) = (b-a)^{-1}$ for $x \in [a,b]$. We deduce that the cumulative density function is equal to:

$$ \mathbf{F}(x) = \begin{cases} 0 & \text{if } x \le a \\ \frac{x-a}{b-a} & \text{if } x \in [a, b] \\ 1 & \text{if } x \ge b \end{cases} $$

<footnote n="10">¹⁰If $p_i = p_j$, we get the Binomial distribution.</footnote>
We have:

$$ \mathbb{E}[X] = \frac{a+b}{2} $$

and:

$$ \mathrm{var}(X) = \frac{(b-a)^2}{12} $$

The standard uniform distribution $U_{[0,1]}$ is obtained with $a = 0$ and $b = 1$.

**Remark 149** The uniform distribution is related to the probability integral transform. We assume that $X \sim \mathbf{F}$ is a continuous random variable. Let $U = \mathbf{F}(X)$ be the integral transform of $X$. Its cumulative distribution function $\mathbf{G}$ is equal to:

$$ \begin{aligned} \mathbf{G}(u) &= \mathrm{Pr}\{U \le u\} \\ &= \mathrm{Pr}\{\mathbf{F}(X) \le u\} \\ &= \mathrm{Pr}\{X \le \mathbf{F}^{-1}(u)\} \\ &= \mathbf{F}(\mathbf{F}^{-1}(u)) \\ &= u \end{aligned} $$

where $\mathbf{G}(0) = 0$ and $\mathbf{G}(1) = 1$. We deduce that $U$ is a standard uniform random variable.

## The gamma distribution

The gamma distribution is a two-parameter family of continuous probability distributions, whose support is $[0, \infty)$. We note $X \sim \mathcal{G}(\alpha, \beta)$ where $\alpha > 0$ and $\beta > 0$. $\alpha$ and $\beta$ are called the shape parameter and the rate parameter. The probability density function is equal to:

$$ f(x) = \frac{\beta^{\alpha} x^{\alpha-1} e^{-\beta x}}{\Gamma(\alpha)} $$

where $\Gamma(\alpha)$ is the gamma function defined as:

$$ \Gamma(\alpha) = \int_{0}^{\infty} t^{\alpha-1} e^{-t} dt $$

The cumulative distribution function is the regularized gamma function:

$$ \mathbf{F}(x) = \frac{\gamma(\alpha, \beta x)}{\Gamma(\alpha)} $$

where $\gamma(\alpha, x)$ is the lower incomplete gamma function defined as:

$$ \gamma(\alpha, x) = \int_{0}^{x} t^{\alpha-1} e^{-t} dt $$

We have $\mathbb{E}[X] = \alpha/\beta$ and $\mathrm{var}(X) = \alpha/\beta^2$. We verify the following properties:

* $\mathcal{G}(1, \beta) \sim \mathcal{E}(\beta)$;

* if $X \sim \mathcal{G}(\alpha, \beta)$, then $cX \sim \mathcal{G}(\alpha, \beta/c)$ when $c > 0$;

* $\sum_{i=1}^{n} \mathcal{G}(\alpha_i, \beta) \sim \mathcal{G}(\sum_{i=1}^{n} \alpha_i, \beta)$.

**Remark 150** The standard gamma distribution corresponds to $\mathcal{G}(\alpha, 1)$ and is denoted by $\mathcal{G}(\alpha)$.
## The beta distribution

The beta distribution is a two-parameter family of continuous probability distributions defined on the interval $[0, 1]$. We note $X \sim \mathcal{B}(\alpha, \beta)$ where $\alpha > 0$ and $\beta > 0$. The probability density function is equal to:

$$f(x) = \frac{x^{\alpha-1} (1-x)^{\beta-1}}{\mathfrak{B}(\alpha, \beta)}$$

where $\mathfrak{B}(\alpha, \beta)$ is the beta function defined as:

$$\mathfrak{B}(\alpha, \beta) = \int_{0}^{1} t^{\alpha-1} (1-t)^{\beta-1} dt = \frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\alpha + \beta)}$$

The cumulative distribution function is the regularized incomplete beta function:

$$\mathbf{F}(x) = \mathcal{IB}(x; \alpha, \beta) = \frac{\mathfrak{B}(x; \alpha, \beta)}{\mathfrak{B}(\alpha, \beta)}$$

where $\mathfrak{B}(x; \alpha, \beta)$ is the incomplete beta function defined as:

$$\mathfrak{B}(x; \alpha, \beta) = \int_{0}^{x} t^{\alpha-1} (1-t)^{\beta-1} dt$$

We have $\mathbb{E}[X] = \alpha / (\alpha + \beta)$ and:

$$\mathrm{var}(X) = \frac{\alpha\beta}{(\alpha + \beta)^2 (\alpha + \beta + 1)}$$

## The exponential distribution

$X$ is an exponential random variable $\mathcal{E}(\lambda)$ if the density function is $f(x) = \lambda e^{-\lambda x}$ for $x \ge 0$. We deduce that $\mathbf{F}(x) = 1 - e^{-\lambda x}$. We have $\mathbb{E}[X] = 1/\lambda$ and $\mathrm{var}(X) = 1/\lambda^2$. More generally, we can show that $\mathbb{E}[X^n] = n!/\lambda^n$. This distribution verifies the lack of memory property:

$$\mathrm{Pr}\{X \ge s + t | X \ge s\} = \mathrm{Pr}\{X \ge t\}$$

for all $s \ge 0$ and $t \ge 0$.

## The normal distribution

Let $\mathbb{C}$ be a correlation matrix. We consider the standardized Gaussian random vector $X \sim \mathcal{N}(\mathbf{0}, \mathbb{C})$ of dimension $n$. We note $\phi_n(x; \mathbb{C})$ the associated density function defined as:

$$\phi_n(x; \mathbb{C}) = (2\pi)^{-n/2} |\mathbb{C}|^{-1/2} \exp \left( -\frac{1}{2} x^\top \mathbb{C}^{-1} x \right)$$

We deduce that the expression of cumulative distribution function is:

$$\Phi_n(x; \mathbb{C}) = \int_{-\infty}^{x_1} \cdots \int_{-\infty}^{x_2} \phi_n(u; \mathbb{C}) du$$

By construction, we have $\mathbb{E}[X] = \mathbf{0}$ and $\mathrm{cov}(x) = \mathbb{C}$. In the bivariate case, we use the notations $\phi_2(x_1, x_2; \rho) = \phi_2(x; \mathbb{C})$ and $\Phi_2(x_1, x_2; \rho) = \Phi_2(x; \mathbb{C})$ where $\rho = \mathbb{C}_{1,2}$ is the correlation between
the components $X_1$ and $X_2$. In the univariate case, we also consider the alternative notations $\phi(x) = \phi_1(x; 1)$ and $\Phi(x) = \Phi_1(x; 1)$. The density function reduces then to:

$$ \phi(x) = \frac{1}{\sqrt{2\pi}} \exp \left( -\frac{1}{2} x^2 \right) $$

Concerning the moments, we have $\mu(X) = 0$, $\sigma(X) = 1$, $\gamma_1(X) = 0$ and $\gamma_2(X) = 0$.

Adding a mean vector $\mu$ and a covariance matrix $\Sigma$ is equivalent to apply the linear transformation to $X$:

$$ Y = \mu + \sigma X $$

where $\sigma = \text{diag}^{1/2}(\Sigma)$.

## The log-normal distribution

Let $Z \sim N(\mu, \sigma^2)$ be a normal-distributed random variable. $X = e^Z$ is a log-normal random variable and we note $X \sim LN(\mu, \sigma^2)$. The probability distribution function is equal to:

$$ f(x) = \frac{1}{x\sigma\sqrt{2\pi}} e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2} $$

whereas the cumulative distribution function has the following expression:

$$ F(x) = \Phi \left( \frac{\ln x - \mu}{\sigma} \right) $$

We have:

$$ \mathbb{E}[X] = e^{\mu + \frac{1}{2}\sigma^2} $$

and:

$$ \text{var}(X) = e^{2\mu + \sigma^2} (e^{\sigma^2} - 1) $$

## The inverse Gaussian distribution

The inverse Gaussian distribution is a two-parameter family of continuous probability distribution, whose support is $[0, \infty)$. We note $X \sim IG(\mu, \lambda)$ where $\mu > 0$ and $\lambda > 0$, $\mu$ and $\lambda$ are called the mean parameter and the shape parameter. The probability density function is equal to:

$$ f(x) = \sqrt{\frac{\lambda}{2\pi x^3}} \exp \left( -\frac{\lambda}{2\mu^2 x} (x - \mu)^2 \right) $$

while the cumulative density function has the following expression:

$$ F(x) = \Phi \left( \sqrt{\frac{\lambda}{x}} \left( \frac{x-\mu}{\mu} \right) \right) + \exp \left( 2\frac{\lambda}{\mu} \right) \Phi \left( -\sqrt{\frac{\lambda}{x}} \left( \frac{x+\mu}{\mu} \right) \right) $$

We have $\mathbb{E}[X] = \mu$ and $\text{var}(X) = \mu^3/\lambda$. We verify the following properties:

* If $X \sim IG(\mu, \lambda)$, then $cX \sim IG(c\mu, c\lambda)$ when $c > 0$;

* If $X_i \sim IG(\mu, \lambda)$, then $\sum_{i=1}^n X_i \sim IG(n\mu, n^2\lambda)$;

* If $X \sim IG(\mu, \lambda)$, then $\mathbb{E}[X^{-k}] = \mu^{-(2k+1)}\mathbb{E}(X^{k+1})$ for $k \in \mathbb{N}^*$.
## The Pareto distribution

The Pareto distribution is denoted by $P(\alpha, x_-)$. We have:

$$f(x) = \frac{\alpha}{x} \left(\frac{x}{x_-}\right)^{-\alpha}$$

and:

$$F(x) = 1 - \left(\frac{x}{x_-}\right)^{-\alpha}$$

where $x \ge x_-$, $\alpha > 0$ and $x_- > 0$. Concerning the first two moments, we obtain:

$$\mathbb{E}[X] = \frac{\alpha x_-}{\alpha - 1}$$

if $\alpha > 1$ and:

$$\mathrm{var}(X) = \frac{\alpha x_-^2}{(\alpha - 1)^2 (\alpha - 2)}$$

if $\alpha > 2$.

**Remark 151** The Pareto distribution: can be parameterized as follows;

$$F(x) = 1 - \left(\frac{\theta + x}{\theta}\right)^{-\alpha}$$

where $x \ge 0$, $\alpha > 0$ and $\theta > 0$. In this case, it is denoted by $P(\alpha, \theta)$.

## The generalized extreme value distribution

The generalized extreme value distribution is denoted by $GEV(\mu, \sigma, \xi)$. We have:

$$f(x) = \frac{1}{\sigma} \left(1 + \xi \left(\frac{x - \mu}{\sigma}\right)\right)^{-(1+1/\xi)} \exp\left(-\left(1 + \xi \left(\frac{x - \mu}{\sigma}\right)\right)^{-1/\xi}\right)$$

and:

$$F(x) = \exp\left(-\left(1 + \xi \left(\frac{x - \mu}{\sigma}\right)\right)^{-1/\xi}\right)$$

where $x > \mu - \sigma/\xi$, $\sigma > 0$ and $\xi > 0$. Concerning the first two moments, we obtain:

$$\mathbb{E}[X] = \mu + \frac{\sigma}{\xi} (\Gamma(1 - \xi) - 1)$$

if $\xi < 1$ and:

$$\mathrm{var}(X) = \frac{\sigma^2}{\xi^2} \left( \Gamma(1 - 2\xi) - \Gamma^2(1 - \xi) \right)$$

if $\xi < 1/2$.
## The generalized Pareto distribution

The generalized Pareto distribution is denoted by $GPD(\sigma, \xi)$. We have:

$$f(x) = \frac{1}{\sigma} \left(1 + \frac{\xi x}{\sigma}\right)^{-1/\xi - 1}$$

and:

$$F(x) = 1 - \left(1 + \frac{\xi x}{\sigma}\right)^{-1/\xi}$$

where $x \ge 0$, $\sigma > 0$ and $\xi > 0$. Concerning the first two moments, we obtain:

$$\mathbb{E}[X] = \frac{\sigma}{1 - \xi}$$

if $\xi < 1$ and:

$$\mathrm{var}(X) = \frac{\sigma^2}{(1 - \xi)^2 (1 - 2\xi)}$$

if $\xi < 1/2$.

## The Bates distribution

Let $U_1, \dots, U_n$ be a sequence of iid uniform random variables $U_{[0,1]}$. The mean $\bar{U}_n$ follows the Bates distribution:

$$X := \bar{U}_n = \frac{1}{n} \sum_{k=1}^{n} U_k \sim \text{Bates}(n)$$

The support is $[0, 1]$. The probability density function is equal to$^{11}$:

$$f(x) = \frac{n}{2(n-1)!} \sum_{k=0}^{n} (-1)^k \binom{n}{k} (nx-k)^{n-1} \operatorname{sgn}(nx-k)$$

whereas the cumulative distribution function has the following expression:

$$F(x) = \frac{1}{n!} \sum_{k=0}^{n} 1 \{nx > k\} (-1)^k \binom{n}{k} (nx-k)^n$$

We have:

$$\mathbb{E}[X] = \frac{1}{2}$$

and

$$\mathrm{var}(X) = \frac{1}{12n}$$

<footnote>$^{11}$If $k = nx$, $\operatorname{sgn}(nx - k) = 0$.</footnote>
## Conditional probability distribution in the Gaussian case

Let us consider a Gaussian random vector defined as follows:

$$ \begin{pmatrix} X \\ Y \end{pmatrix} \sim \mathcal{N} \left( \begin{pmatrix} \mu_x \\ \mu_y \end{pmatrix}, \begin{pmatrix} \Sigma_{x,x} & \Sigma_{x,y} \\ \Sigma_{y,x} & \Sigma_{y,y} \end{pmatrix} \right) $$

The conditional probability distribution of $Y$ given $X = x$ is a multivariate normal distribution. We have:

$$ \mu_{y|x} = \mathbb{E}[Y | X = x] = \mu_y + \Sigma_{y,x} \Sigma_{x,x}^{-1} (x - \mu_x) $$

and:

$$ \Sigma_{y,y|x} = \sigma^2 [Y | X = x] = \Sigma_{y,y} - \Sigma_{y,x} \Sigma_{x,x}^{-1} \Sigma_{x,y} $$

We deduce that:

$$ Y = \mu_y + \Sigma_{y,x} \Sigma_{x,x}^{-1} (x - \mu_x) + u $$

where $u$ is a centered Gaussian random variable with variance $\sigma^2 = \Sigma_{y,y|x}$. It follows that:

$$ Y = \underbrace{(\mu_y - \Sigma_{y,x} \Sigma_{x,x}^{-1} \mu_x)}_{\beta_0} + \underbrace{\Sigma_{y,x} \Sigma_{x,x}^{-1} x}_{\beta^\top} + u $$

We recognize the linear regression of $Y$ on a constant and a set of exogenous variables $X$:

$$ Y = \beta_0 + \beta^\top X + u $$

Moreover, we have:

$$ \mathfrak{R}^2 = 1 - \frac{\sigma^2}{\Sigma_{y,y}} = \frac{\Sigma_{y,x} \Sigma_{x,x}^{-1} \Sigma_{x,y}}{\Sigma_{y,y}} $$

## Change of variables

Let $X$ be a random variable whose probability density function is $f(x)$. We consider the change of the variable $Y = \varphi(X)$. If the function $\varphi$ is monotone, then the probability density function $g(y)$ of $Y$ is equal to:

$$ g(y) = f(x) \left| \frac{dx}{dy} \right| $$

In the multivariate case, let $(X_1, \dots, X_n)$ be the random vector with density function $f(x_1, \dots, x_n)$. If the function $\varphi$ is bijective, we can show that the probability density function of $(Y_1, \dots, Y_n) = \varphi(X_1, \dots, X_n)$ is equal to:

$$ g(y_1, \dots, y_n) = f(x_1, \dots, x_n) \left| \frac{1}{\det J_{\varphi}} \right| $$

where $J_{\varphi}$ is the Jacobian associated with the variable change.
Ratio distribution

Given two independent random variables X and Y, the distribution of the random variable Z = Y/X is a ratio distribution. We have:

$$
\begin{align*}
\mathbf{F}_z(z) &= \operatorname{Pr}\{Z \le z\} \\
&= \operatorname{Pr}\left\{\frac{X}{Y} \le z\right\} \\
&= \operatorname{Pr}\{X \ge zY, Y \le 0\} + \operatorname{Pr}\{X \le zY, Y > 0\} \\
&= \int_{-\infty}^{0} \int_{zy}^{\infty} f_{x,y}(x,y) \, dx \, dy + \int_{0}^{\infty} \int_{-\infty}^{zy} f_{x,y}(x,y) \, dx \, dy \\
&= \int_{-\infty}^{0} \int_{zy}^{\infty} f_x(x) f_y(y) \, dx \, dy + \int_{0}^{\infty} \int_{-\infty}^{zy} f_x(x) f_y(y) \, dx \, dy \\
&= \int_{-\infty}^{0} \left( \int_{zy}^{\infty} f_x(x) \, dx \right) f_y(y) \, dy + \int_{0}^{\infty} \left( \int_{-\infty}^{zy} f_x(x) \, dx \right) f_y(y) \, dy \\
&= \int_{-\infty}^{0} (1 - \mathbf{F}_x(zy)) f_y(y) \, dy + \int_{0}^{\infty} \mathbf{F}_x(zy) f_y(y) \, dy
\end{align*}
$$

and:

$$
\begin{align*}
f_z(z) &= \partial_z \operatorname{Pr}\{Z \le z\} \\
&= -\int_{-\infty}^{0} y f_x(zy) f_y(y) \, dy + \int_{0}^{\infty} y f_x(zy) f_y(y) \, dy \\
&= \int_{-\infty}^{0} |y| f_x(zy) f_y(y) \, dy + \int_{0}^{\infty} |y| f_x(zy) f_y(y) \, dy \\
&= \int_{-\infty}^{\infty} |y| f_x(zy) f_y(y) \, dy
\end{align*}
$$

In the case where $X \sim \mathcal{N}(\mu_x, \sigma_x^2)$ and $Y \sim \mathcal{N}(\mu_y, \sigma_y^2)$, we note $Z \sim \mathcal{NRD}(\mu_x, \sigma_x^2, \mu_y, \sigma_y^2)$, and Hinkley (1969) showed that:

$$
\mathbf{F}_z(z) = \Phi_2 \left( -\frac{\mu_y z - \mu_x}{\sigma_x \sigma_y a(z)}, -\frac{\mu_y}{\sigma_y}, \rho_z \right) + \Phi_2 \left( \frac{\mu_y z - \mu_x}{\sigma_x \sigma_y a(z)}, \frac{\mu_y}{\sigma_y}, \rho_z \right)
$$

and:

$$
f_z(z) = \frac{b(z)}{\sigma_x \sigma_y \sqrt{2\pi a^3(z)}} \left( 2\Phi\left(\frac{b(z)}{a(z)}\right) - 1 \right) \exp\left(\frac{b^2(z) - ca^2(z)}{2a^2(z)}\right) + \frac{1}{\sigma_x \sigma_y a^2(z)} \pi \exp\left(-\frac{c}{2}\right)
$$

where $a(z) = \sqrt{\frac{1}{\sigma_x^2}z^2 + \frac{1}{\sigma_y^2}}$, $b(z) = \frac{\mu_x}{\sigma_x^2}z + \frac{\mu_y}{\sigma_y^2}$, $c = \frac{\mu_x^2}{\sigma_x^2} + \frac{\mu_y^2}{\sigma_y^2}$ and $\rho_z = \frac{z}{\sigma_x a(z)}$. Depending on the parameter values, Z has one or two modes. The moments of Z do not exist. In the case where $\mu_x = \mu_y = 0$, $b(z) = 0$ and $c = 0$, which implies that Z follows a Cauchy distribution:

$$
f_z(z) = \frac{\sigma_x \sigma_y}{\pi (\sigma_y^2 z^2 + \sigma_x^2)}
$$

The previous result can be easily extended to the correlated case, because $\mathbf{F}_z(z)$ has exactly the same expression with $a(z) = \sqrt{\frac{1}{\sigma_x^2}z^2 - 2\frac{\rho_{x,y}}{\sigma_x\sigma_y}z + \frac{1}{\sigma_y^2}}$ and $\rho_z = \frac{\sigma_y z - \rho_{x,y}\sigma_x}{\sigma_x\sigma_y a(z)}$.
## A.2.2 Copula functions

The concept of copula has been introduced by Sklar in 1959. During a long time, only a small number of people have used copula functions, more in the field of mathematics than this of statistics. The publication of Genest and MacKay (1986b) in the American Statistician marks a breakdown and opens areas of study in empirical modeling, statistics and econometrics. In what follows, we intensively use the materials developed in the books of Joe (1997), Nelsen (2006) and Roncalli (2020a).

### Definition and main properties

Nelsen (2006) defines a bi-dimensional copula (or a 2-copula) as a function C which satisfies the following properties:

1. Dom C = [0, 1] × [0, 1];

2. C(0, u) = C(u, 0) = 0 and C(1, u) = C(u, 1) = u for all u in [0, 1];

3. C is 2-increasing:

$$C(v_1, v_2) - C(v_1, u_2) - C(u_1, v_2) + C(u_1, u_2) \geq 0$$

for all $(u_1, u_2) \in [0, 1]^2$, $(v_1, v_2) \in [0, 1]^2$ such that $0 \le u_1 \le v_1 \le 1$ and $0 \le u_2 \le v_2 \le 1$.

This definition means that C is a cumulative distribution function with uniform margins:

$$C(u_1, u_2) = \Pr\{U_1 \le u_1, U_2 \le u_2\}$$

where $U_1$ and $U_2$ are two uniform random variables.

Let us consider the function $C^{\perp}(u_1, u_2) = u_1 u_2$. We have $C^{\perp}(0, u) = C^{\perp}(u, 0) = 0$ and $C^{\perp}(1, u) = C^{\perp}(u, 1) = u$. Since we have $v_2 - u_2 \ge 0$ and $v_1 \ge u_1$, it follows that $v_1(v_2 - u_2) \ge u_1(v_2 - u_2)$ and $v_1 v_2 + u_1 u_2 - u_1 v_2 - v_1 u_2 \ge 0$. We deduce that $C^{\perp}$ is a copula function. It is called the product copula.

Let $F_1$ and $F_2$ be any two univariate distributions. It is obvious that $F(x_1, x_2) = C(F_1(x_1), F_2(x_2))$ is a probability distribution with margins $F_1$ and $F_2$. Indeed, $u_i = F_i(x_i)$ defines a uniform transformation $(u_i \in [0, 1])$. Moreover, we verify that $C(F_1(x_1), F_2(\infty)) = C(F_1(x_1), 1) = F_1(x_1)$. Copulas are then a powerful tool to build a multivariate probability distribution when the margins are given. Conversely, Sklar (1959) proves that any bivariate distribution $F$ admits such a representation:

$$F(x_1, x_2) = C(F_1(x_1), F_2(x_2)) \qquad (A.11)$$

and that the copula C is unique provided the margins are continuous. This result is important, because we can associate to each bivariate distribution a copula function. We then obtain a canonical representation of a bivariate probability distribution: on one side, we have the margins or the univariate directions $F_1$ and $F_2$; on the other side, we have the copula C that links these margins and gives the dependence between the unidimensional directions.

If the joint distribution function $F(x_1, x_2)$ is absolutely continuous, we obtain:

$$ \begin{aligned} f(x_1, x_2) &= \partial_{1,2} F(x_1, x_2) \\ &= \partial_{1,2} C(F_1(x_1), F_2(x_2)) \\ &= c(F_1(x_1), F_2(x_2)) \cdot f_1(x_1) \cdot f_2(x_2) \end{aligned} \qquad (A.12) $$
where $f(x_1, x_2)$ is the joint probability density function, $f_1$ and $f_2$ are the marginal densities and $c$ is the copula density:

$$c(u_1, u_2) = \partial_{1,2} \mathbf{C}(u_1, u_2)$$

We notice that the condition $\mathbf{C}(v_1, v_2) - \mathbf{C}(v_1, u_2) - \mathbf{C}(u_1, v_2) + \mathbf{C}(u_1, u_2) \ge 0$ is then equivalent to $\partial_{1,2} \mathbf{C}(u_1, u_2) \ge 0$ when the copula density exists. From Equation (A.12), we deduce that:

$$c(u_1, u_2) = \frac{f(\mathbf{F}_1^{-1}(u_1), \mathbf{F}_2^{-1}(u_2))}{f_1(\mathbf{F}_1^{-1}(u_1)) \cdot f_2(\mathbf{F}_2^{-1}(u_2))} \quad (A.13)$$

We obtain a second canonical representation based on density functions. For some copulas, there is no explicit analytical formula. This is the case of the Normal copula, which is equal to $\mathbf{C}(u_1, u_2; \rho) = \Phi(\Phi^{-1}(u_1), \Phi^{-1}(u_2); \rho)$. Using Equation (A.13), we can however characterize its density function:

$$c(u_1, u_2; \rho) = \frac{2\pi (1 - \rho^2)^{-1/2} \exp \left( -\frac{1}{2(1-\rho^2)} (x_1^2 + x_2^2 - 2\rho x_1 x_2) \right)}{(2\pi)^{-1/2} \exp \left( -\frac{1}{2} x_1^2 \right) \cdot (2\pi)^{-1/2} \exp \left( -\frac{1}{2} x_2^2 \right)} \\ = \frac{1}{\sqrt{1-\rho^2}} \exp \left( -\frac{1}{2} \frac{(x_1^2 + x_2^2 - 2\rho x_1 x_2)}{(1-\rho^2)} + \frac{1}{2} (x_1^2 + x_2^2) \right)$$

where $x_1 = \mathbf{F}_1^{-1}(u_1)$ and $x_2 = \mathbf{F}_2^{-1}(u_2)$. It is then easy to generate bivariate non-normal distributions.

**Example 68** In Figure A.1, we have built a bivariate probability distribution by considering that the margins are a Student's distribution and a beta distribution. The copula function corresponds to the Normal copula such that its Kendall's tau is equal to 60%.

## Fréchet classes and concordance ordering

The goal of Fréchet classes is to study the structure of the class of distributions with given margins. Let us first consider the bivariate case. The distribution function $\mathbf{F}$ belongs to the Fréchet class $(\mathbf{F}_1, \mathbf{F}_2)$ and we note $\mathbf{F} \in \mathcal{F}(\mathbf{F}_1, \mathbf{F}_2)$ if an only if the margins of $\mathbf{F}$ are $\mathbf{F}_1$ and $\mathbf{F}_2$, meaning that $\mathbf{F}(x_1, \infty) = \mathbf{F}_1(x_1)$ and $\mathbf{F}(\infty, x_2) = \mathbf{F}_2(x_2)$. Characterizing the Fréchet class $\mathcal{F}(\mathbf{F}_1, \mathbf{F}_2)$ is then equivalent to find the set $\mathcal{C}$ of copula functions:

$$\mathcal{F}(\mathbf{F}_1, \mathbf{F}_2) = \{ \mathbf{F} : \mathbf{F}(x_1, x_2) = \mathbf{C}(\mathbf{F}_1(x_1), \mathbf{F}_2(x_2)), \mathbf{C} \in \mathcal{C} \}$$

Therefore this problem does not depend on the margins $\mathbf{F}_1$ and $\mathbf{F}_2$. We can show that the extremal distribution functions $\mathbf{F}^-$ and $\mathbf{F}^+$ of the Fréchet class $\mathcal{F}(\mathbf{F}_1, \mathbf{F}_2)$ are $\mathbf{F}^-(x_1, x_2) = \max(\mathbf{F}_1(x_1) + \mathbf{F}_2(x_2) - 1, 0)$ and $\mathbf{F}^+(x_1, x_2) = \min(\mathbf{F}_1(x_1), \mathbf{F}_2(x_2))$. $\mathbf{F}^-$ and $\mathbf{F}^+$ are called the Fréchet lower and upper bounds. We deduce that the corresponding copula functions are:

$$\mathbf{C}^{-}(u_1, u_2) = \max(u_1 + u_2 - 1, 0)$$

and:

$$\mathbf{C}^{+}(u_1, u_2) = \min(u_1, u_2)$$
<figure><img src="image_533.png" /><figcaption>Figure A.1: Example of a bivariate probability distribution with given margins</figcaption></figure>

The extension of bivariate copulas to multivariate copulas is straightforward. Thus, the canonical decomposition of a multivariate distribution function is:

$$F(x_1, \dots, x_n) = C(F_1(x_1), \dots, F_n(x_n))$$

We note $C_E$ the sub-copula of $C$ such that arguments that are not in the set $\mathcal{E}$ are equal to 1. For instance, with a dimension of 4, we have $C_{12}(u, v) = C(u, v, 1, 1)$ and $C_{124}(u, v, w) = C(u, v, 1, w)$. Let us consider the 2-copulas $C_1$ and $C_2$. It seems logical to build a copula of higher dimension with copulas of lower dimensions. In fact, the function $C_1(u_1, C_2(u_2, u_3))$ is not a copula in most cases. In the multivariate case, we define:

$$C^{-}(u_1, \dots, u_n) = \max \left( \sum_{i=1}^{n} u_i - n + 1, 0 \right)$$

and:

$$C^{+}(u_1, \dots, u_n) = \min(u_1, \dots, u_n)$$

We can show that $C^+$ is a copula, but $C^-$ does not belong to the set $\mathcal{C}$. Nevertheless, $C^-$ is the best-possible bound, meaning that for all $(u_1, \dots, u_n) \in [0, 1]^n$, there is a copula that coincide with $C^-$ (Nelsen, 2006). This implies that $F(F_1, \dots, F_n)$ has a minimal distribution function if and only if $\max(\sum_{i=1}^n F_i(x_i) - n + 1, 0)$ is a probability distribution.

We now introduce a stochastic ordering on copulas. Let $C_1$ and $C_2$ be two copula functions. We say that the copula $C_1$ is smaller than the copula $C_2$ and we note $C_1 \prec C_2$ if we verify that $C_1(u_1, u_2) \le C_2(u_1, u_2)$ for all $(u_1, u_2) \in [0, 1]^2$. This stochastic ordering is called the concordance ordering and may be viewed as the first order of the stochastic dominance on probability
distributions. Using the previous results on Fréchet classes, we deduce that:

$$C^{-} \prec C \prec C^{+}$$

for all $C \in \mathcal{C}$. It follows that $C^{-} \prec C^{\perp} \prec C^{+}$. A copula $C$ has a positive quadrant dependence (PQD) if it satisfies the inequality $C^{\perp} \prec C \prec C^{+}$. In a similar way, $C$ has a negative quadrant dependence (NQD) if it satisfies the inequality $C^{-} \prec C \prec C^{\perp}$. As it is a partial ordering, there exist copula functions $C$ such that $C \nprec C^{\perp}$ and $C \nprec C^{\perp}$. A copula function may then have a dependence structure that is neither positive or negative. In Figure A.2, we report the cumulative distribution function (above panel) and its contour lines (below panel) of the three copula functions $C^{-}$, $C^{\perp}$ and $C^{+}$, which plays an important role to understand the dependence between unidimensional margins.

<figure><img src="image_534.png" /><figcaption>Figure A.2: The three copula functions $C^{-}$, $C^{\perp}$ and $C^{+}$</figcaption></figure>

Let $C_{\theta}(u_1, u_2) = C(u_1, u_2; \theta)$ be a family of copula functions that depends on the parameter $\theta$. The copula family $\{C_{\theta}\}$ is totally ordered if, for all $\theta_2 \ge \theta_1$, $C_{\theta_2} \succ C_{\theta_1}$ (positively ordered) or $C_{\theta_2} \prec C_{\theta_1}$ (negatively ordered). For instance, the Frank copula defined by:

$$C(u_1, u_2; \theta) = -\frac{1}{\theta} \ln \left( 1 + \frac{(e^{-\theta u_1} - 1)(e^{-\theta u_2} - 1)}{e^{-\theta} - 1} \right)$$

where $\theta \in \mathbb{R}$ is a positively ordered family. An illustration is showed in Figure A.3). For a given value $\alpha \in [0, 1]$, we also verify that the level curves of $C$ are in the triangle defined as follows:

$$\{(u_1, u_2) : \max(u_1 + u_2 - 1, 0) \le \alpha, \min(u_1, u_2) \ge \alpha\}$$
<figure><img src="image_535.png" /><figcaption>Figure A.3: Concordance ordering of the Frank copula</figcaption></figure>

## Copula function and random vector

Let $X = (X_1, X_2)$ be a random vector with distribution $\mathbf{F}$. We define the copula of $(X_1, X_2)$ by the copula of $\mathbf{F}$:

$$ \mathbf{F}(x_1, x_2) = \mathbf{C}\langle X_1, X_2 \rangle (\mathbf{F}_1(x_1), \mathbf{F}_2(x_2)) $$

In what follows, we give the main results on the dependence of the random vector $X$ found in Nelsen (2006). We first consider the probabilistic interpretation of the three copula functions $\mathbf{C}^{-}$, $\mathbf{C}^{\perp}$ and $\mathbf{C}^{+}$:

* $X_1$ and $X_2$ are countermonotonic — or $\mathbf{C}\langle X_1, X_2 \rangle = \mathbf{C}^{-}$ — if there exists a random variable $X$ such that $X_1 = f_1(X)$ and $X_2 = f_2(X)$ where $f_1$ and $f_2$ are respectively decreasing and increasing functions<fnref n="12" />;

* $X_1$ and $X_2$ are independent if the dependence function is the product copula $\mathbf{C}^{\perp}$;

* $X_1$ are $X_2$ are comonotonic — or $\mathbf{C}\langle X_1, X_2 \rangle = \mathbf{C}^{+}$ — if there exists a random variable $X$ such that $X_1 = f_1(X)$ and $X_2 = f_2(X)$ where $f_1$ and $f_2$ are both increasing functions<fnref n="13" />.

Let us consider a uniform random vector $(U_1, U_2)$. We have $U_2 = 1 - U_1$ when $\mathbf{C}\langle X_1, X_2 \rangle = \mathbf{C}^{-}$ and $U_2 = U_1$ when $\mathbf{C}\langle X_1, X_2 \rangle = \mathbf{C}^{+}$. In the case of a standardized Gaussian random vector, we obtain $X_2 = -X_1$ when $\mathbf{C}\langle X_1, X_2 \rangle = \mathbf{C}^{-}$ and $X_2 = X_1$ when $\mathbf{C}\langle X_1, X_2 \rangle = \mathbf{C}^{+}$. If the marginals are log-normal, it follows that $X_2 = X_1^{-1}$ when $\mathbf{C}\langle X_1, X_2 \rangle = \mathbf{C}^{-}$ and $X_2 = X_1$ when $\mathbf{C}\langle X_1, X_2 \rangle = \mathbf{C}^{+}$. For these three examples, we verify that $X_2$ is a decreasing (resp. increasing) function of $X_1$ if

<footnote n="12">¹²We also have $X_2 = f(X_1)$ where $f = f_2 \circ f_1^{-1}$ is a decreasing function.</footnote>

<footnote n="13">¹³In this case, $X_2 = f(X_1)$ where $f = f_2 \circ f_1^{-1}$ is an increasing function.</footnote>
the copula function $C(X_1, X_2)$ is $C^-$ (resp. $C^+$). The concepts of counter- and comonotonicity concepts generalize the cases where the linear correlation of a Gaussian vector is equal to -1 or +1. Indeed, $C^-$ and $C^+$ define respectively perfect negative and positive dependence.

Let $(X_1, X_2)$ be a random vectors, whose copula is $C(X_1, X_2)$. If $h_1$ and $h_2$ are two increasing functions on $\text{Im } X_1$ and $\text{Im } X_2$, then we have:

$$C(h_1(X_1), h_2(X_2)) = C(X_1, X_2)$$

This means that copula functions are invariant under strictly increasing transformations of the random variables. To prove this theorem, we note $\mathbf{F}$ and $\mathbf{G}$ the probability distributions of the random vectors $(X_1, X_2)$ and $(Y_1, Y_2) = (h_1(X_1), h_2(X_2))$. The margins of $\mathbf{G}$ are:

$$ 
\begin{aligned}
\mathbf{G}_1 (y_1) &= \Pr \{Y_1 \le y_1\} \\
&= \Pr \{h_1 (X_1) \le y_1\} \\
&= \Pr \{X_1 \le h_1^{-1} (y_1)\} \quad (\text{because } h_1 \text{ is strictly increasing}) \\
&= \mathbf{F}_1 (h_1^{-1} (y_1))
\end{aligned}
$$

and $\mathbf{G}_2 (y_2) = \mathbf{F}_2 (h_2^{-1} (y_2))$. We deduce that $\mathbf{G}_1^{-1} (u_1) = h_1 (\mathbf{F}_1^{-1} (u_1))$ and $\mathbf{G}_2^{-1} (u_2) = h_2 (\mathbf{F}_2^{-1} (u_2))$. By definition, we have:

$$C(Y_1, Y_2) (u_1, u_2) = \mathbf{G} (\mathbf{G}_1^{-1} (u_1), \mathbf{G}_2^{-1} (u_2))$$

Moreover, it follows that:

$$ 
\begin{aligned}
\mathbf{G} (\mathbf{G}_1^{-1} (u_1), \mathbf{G}_2^{-1} (u_2)) &= \Pr \{Y_1 \le \mathbf{G}_1^{-1} (u_1), Y_2 \le \mathbf{G}_2^{-1} (u_2)\} \\
&= \Pr \{h_1 (X_1) \le \mathbf{G}_1^{-1} (u_1), h_2 (X_2) \le \mathbf{G}_2^{-1} (u_2)\} \\
&= \Pr \{X_1 \le h_1^{-1} (\mathbf{G}_1^{-1} (u_1)), X_2 \le h_2^{-1} (\mathbf{G}_2^{-1} (u_2))\} \\
&= \Pr \{X_1 \le \mathbf{F}_1^{-1} (u_1), X_2 \le \mathbf{F}_2^{-1} (u_2)\} \\
&= \mathbf{F} (\mathbf{F}_1^{-1} (u_1), \mathbf{F}_2^{-1} (u_2))
\end{aligned}
$$

Because we have $C(X_1, X_2) (u_1, u_2) = \mathbf{F} (\mathbf{F}_1^{-1} (u_1), \mathbf{F}_2^{-1} (u_2))$, we deduce that $C(Y_1, Y_2) = C(X_1, X_2)$.

We can interpret the copula function $C(X_1, X_2)$ as a standardization of the joint distribution after eliminating the effects of margins. Indeed, it is a comprehensive statistic of the dependence function between $X_1$ and $X_2$. Therefore, a non-comprehensive statistic will be a dependence measure if it can be expressed using $C(X_1, X_2)$. Following Nelsen (2006), a numeric measure $m$ of association between $X_1$ and $X_2$ is a measure of concordance if it satisfies the following properties:

1. $-1 = m(X, -X) \le m(C) \le m(X, X) = 1$;

2. $m(C^\perp) = 0$;

3. $m(-X_1, X_2) = m(X_1, -X_2) = -m(X_1, X_2)$;

4. if $C_1 \prec C_2$, then $m(C_1) \le m(C_2)$;

Using this last property, we have: $C \prec C^\perp \Rightarrow m(C) < 0$ and $C \succ C^\perp \Rightarrow m(C) > 0$. The concordance measure can then be viewed as a generalization of the linear correlation when the dependence function is not normal. Indeed, a positive quadrant dependence copula will have a positive concordance measure whereas a negative quadrant dependence copula will have a negative
concordance measure. Moreover, the bounds -1 and +1 are reached when the copula function is countermonotonic and comonotonic.

Among the several concordance measures, we find Kendall's tau and Spearman's rho, which play an important role in non-parametric statistics. Let us consider a sample of n observations {($x_1, y_1$), ..., ($x_n, y_n$)} of the random vector (X, Y). Kendall's tau is the probability of concordance — ($X_i - X_j$) $\cdot$ ($Y_i - Y_j$) > 0 — minus the probability of discordance — ($X_i - X_j$) $\cdot$ ($Y_i - Y_j$) < 0:

$$ \tau = \Pr \{(X_i - X_j) \cdot (Y_i - Y_j) > 0\} - \Pr \{(X_i - X_j) \cdot (Y_i - Y_j) < 0\} $$

Spearman's rho is the linear correlation of the rank statistics ($X_{i:n}, Y_{i:n}$). We can also show that Spearman's rho has the following expression:

$$ \varrho = \frac{\operatorname{cov}(\mathbf{F}_X(X), \mathbf{F}_Y(Y))}{\sigma(\mathbf{F}_X(X)) \cdot \sigma(\mathbf{F}_Y(Y))} $$

Schweizer and Wolff (1981) showed that Kendall's tau and Spearman's rho are concordance measures and have the following expressions:

$$ \tau = 4 \iint_{[0,1]^2} \mathbf{C}(u_1, u_2) \, d\mathbf{C}(u_1, u_2) - 1 $$

$$ \varrho = 12 \iint_{[0,1]^2} u_1 u_2 \, d\mathbf{C}(u_1, u_2) - 3 $$

From a numerical point of view, the following formulas should be preferred (Nelsen, 2006):

$$ \tau = 1 - 4 \iint_{[0,1]^2} \partial_{u_1} \mathbf{C}(u_1, u_2) \partial_{u_2} \mathbf{C}(u_1, u_2) \, du_1 du_2 $$

$$ \varrho = 12 \iint_{[0,1]^2} \mathbf{C}(u_1, u_2) \, du_1 du_2 - 3 $$

For some copulas, we have analytical formulas. For instance, we have:

<table><thead><tr><th>Copula</th><th>ρ</th><th>τ</th></tr></thead><tbody><tr><td>Normal</td><td>6π<sup>-1</sup> arc sin (ρ/2)</td><td>2π<sup>-1</sup> arc sin (ρ)</td></tr><tr><td>Gumbel</td><td>✓</td><td>(θ - 1) / θ</td></tr><tr><td>FGM</td><td>θ/3</td><td>2θ/9</td></tr><tr><td>Frank</td><td>1 - 12θ<sup>-1</sup> (D<sub>1</sub>(θ) - D<sub>2</sub>(θ))</td><td>1 - 4θ<sup>-1</sup> (1 - D<sub>1</sub>(θ))</td></tr></tbody></table>

where $D_k(x)$ is the Debye function. The Gumbel (or Gumbel-Hougaard) copula is equal to:

$$ \mathbf{C}(u_1, u_2; \theta) = \exp \left( - \left[ (-\ln u_1)^\theta + (-\ln u_2)^\theta \right]^{1/\theta} \right) $$

for $\theta \ge 1$, whereas the expression of the Farlie-Gumbel-Morgenstern (or FGM) copula is:

$$ \mathbf{C}(u_1, u_2; \theta) = u_1 u_2 (1 + \theta (1 - u_1) (1 - u_2)) $$

for $-1 \le \theta \le 1$.

For illustration, we report in Figures A.4, A.5 and A.6 the level curves of several density functions built with Normal, Frank and Gumbel copulas. In order to compare them, the parameter of each copula is calibrated such that Kendall's tau is equal to 50%. This means that these 12 distributions functions have the same dependence with respect to Kendall's tau. However, the dependence is different from one figure to another, because their copula function is not the same. This is why Kendall's tau is not an exhaustive statistic of the dependence between two random variables.
<figure><img src="image_536.png" /><figcaption>Figure A.5: Contour lines of bivariate densities (Frank copula)</figcaption></figure>

<figure><img src="image_537.png" /></figure>
<figure><img src="image_538.png" /><figcaption>Figure A.6: Contour lines of bivariate densities (Gumbel copula)</figcaption></figure>

**Remark 152** We can show that the linear correlation (or Pearson's correlation) is not a concordance measure (Roncalli, 2020a, pages 727-729). In particular, the lower and upper bounds $\rho \langle C^- \rangle$ and $\rho \langle C^+ \rangle$ are not necessarily equal to -1 and +1. While the copula function is an invariant measure by increasing transformations, the Pearson's correlation is an invariant measure by increasing linear transformations. This is why the correlation is called a linear dependence measure.

## Parametric copula functions

**Archimedean copulas** Genest and MacKay (1986b) define Archimedean copulas as follows:

$$ C(u_1, u_2) = \begin{cases} \varphi^{-1}(\varphi(u_1) + \varphi(u_2)) & \text{if } \varphi(u_1) + \varphi(u_2) \le \varphi(0) \\ 0 & \text{otherwise} \end{cases} $$

where $\varphi$ is a $C^2$ is a function which satisfies $\varphi(1) = 0$, $\varphi'(u) < 0$ and $\varphi''(u) > 0$ for all $u \in [0, 1]$. $\varphi(u)$ is called the generator of the copula function. If $\varphi(0) = \infty$, the generator is said to be strict. Genest and MacKay (1986a) MacKay, Jock link the construction of Archimedean copulas to the independence of random variables. Indeed, by considering the multiplicative generator $\lambda(u) = \exp(-\varphi(u))$, the authors show that:

$$ C(u_1, u_2) = \lambda^{-1}(\lambda(u_1) \lambda(u_2)) $$

This means that:

$$ \lambda(\Pr\{U_1 \le u_1, U_2 \le u_2\}) = \lambda(\Pr\{U_1 \le u_1\}) \times \lambda(\Pr\{U_2 \le u_2\}) $$

In this case, the random variables $(U_1, U_2)$ become independent when the scale of probabilities has been transformed.
The product copula $\mathbf{C}^{\perp}$ is Archimedean and the associated generator is $\varphi(u) = -\ln u$. Concerning Fréchet copulas, only $\mathbf{C}^{-}$ is Archimedean with $\varphi(u) = 1 - u$. In Table A.1, we provide another examples of Archimedean copulas$^{14}$.

<table><caption>Table A.1: Archimedean copula functions</caption><thead><tr><th>Copula</th><th>&phi; (u)</th><th>C (u<sub>1</sub>, u<sub>2</sub>)</th></tr></thead><tbody><tr><td>Clayton</td><td>u<sup>-&theta;</sup> - 1</td><td>(u<sub>1</sub><sup>-&theta;</sup> + u<sub>2</sub><sup>-&theta;</sup> - 1)<sup>-1/&theta;</sup></td></tr><tr><td>Frank</td><td>-&ln; (e<sup>-&theta;u</sup> - 1) / (e<sup>-&theta;</sup> - 1)</td><td>-1/&theta; &ln; (1 + (e<sup>-&theta;u<sub>1</sub></sup> - 1) / (e<sup>-&theta;u<sub>2</sub></sup> - 1))</td></tr><tr><td>Gumbel</td><td>(-&ln; u)<sup>&theta;</sup></td><td>&exp; (-(u&#771;<sub>1</sub><sup>&theta;</sup> + u&#771;<sub>2</sub><sup>&theta;</sup>)<sup>1/&theta;</sup>)</td></tr><tr><td>Joe</td><td>-&ln; (1 - (1 - u)<sup>&theta;</sup>)</td><td>1 - (u&#771;<sub>1</sub><sup>&theta;</sup> + u&#771;<sub>2</sub><sup>&theta;</sup> - u&#771;<sub>1</sub><sup>&theta;</sup>u&#771;<sub>2</sub><sup>&theta;</sup>)<sup>1/&theta;</sup></td></tr></tbody></table>

Nelsen (2006) showed that if $\varphi(t)$ is a strict generator, then we can build two-parameter Archimedean copulas by considering the following generator:

$$ \varphi_{\alpha, \beta}(t) = (\varphi(t^{\alpha}))^{\beta} $$

where $\alpha > 0$ and $\beta > 1$. For instance, if $\varphi(t) = t^{-1} - 1$, the two-parameter generator is $\varphi_{\alpha, \beta}(t) = (t^{-\alpha} - 1)^{\beta}$. Therefore, the corresponding copula function is defined by:

$$ \mathbf{C}(u_1, u_2) = \left( \left[ (u_1^{-\alpha} - 1)^{\beta} + (u_2^{-\alpha} - 1)^{\beta} \right]^{1/\beta} + 1 \right)^{-1/\alpha} $$

This is a generalization of the Clayton copula, which is obtained when the parameter $\beta$ is equal to 1.

We can build multivariate Archimedean copulas in the following way:

$$ \mathbf{C}(u_1, \dots, u_n) = \varphi^{-1} (\varphi(u_1) + \dots + \varphi(u_n)) $$

However, $\mathbf{C}$ is a copula function if and only if the function $\varphi^{-1}(u)$ is completely monotone (Nelsen, 2006):

$$ (-1)^k \frac{d^k}{du^k} \varphi^{-1}(u) \geq 0 \quad \forall k \geq 1 $$

For instance, the multivariate GumbelCopula function!Gumbel-Hougaard copula is defined by:

$$ \mathbf{C}(u_1, \dots, u_n) = \exp \left( - \left( (-\ln u_1)^\theta + \dots + (-\ln u_n)^\theta \right)^{1/\theta} \right) $$

**Normal copula** The Normal copula is the dependence function of the multivariate normal distribution with a correlation matrix $\rho$:

$$ \mathbf{C}(u_1, \dots, u_n; \rho) = \Phi_n (\Phi^{-1}(u_1), \dots, \Phi^{-1}(u_n); \rho) $$

By using the canonical decomposition of the multivariate density function:

$$ f(x_1, \dots, x_n) = c(\mathbf{F}_1(x_1), \dots, \mathbf{F}_n(x_n)) \prod_{i=1}^{n} f_i(x_i) $$

<footnote n="14">¹⁴We use the notations $\bar{u} = 1 - u$ and $\tilde{u} = -\ln u$.</footnote>
we deduce that the probability density function of the Normal copula is:

$$c(u_1, \dots, u_n; \rho) = \frac{1}{|\rho|^{\frac{1}{2}}} \exp \left( -\frac{1}{2} x^{\top} (\rho^{-1} - I_n) x \right)$$

where $x_i = \Phi^{-1}(u_i)$. In the bivariate case, we obtain<fnref n="15" />:

$$c(u_1, u_2; \rho) = \frac{1}{\sqrt{1 - \rho^2}} \exp \left( - \frac{x_1^2 + x_2^2 - 2\rho x_1 x_2}{2(1 - \rho^2)} + \frac{x_1^2 + x_2^2}{2} \right)$$

It follows that the expression of the bivariate Normal copula function is also equal to:

$$C(u_1, u_2; \rho) = \int_{-\infty}^{\Phi^{-1}(u_1)} \int_{-\infty}^{\Phi^{-1}(u_2)} \phi_2(x_1, x_2; \rho) \, dx_1 \, dx_2$$

where $\phi_2(x_1, x_2; \rho)$ is the bivariate normal density:

$$\phi_2(x_1, x_2; \rho) = \frac{1}{2\pi\sqrt{1-\rho^2}} \exp\left(-\frac{x_1^2 + x_2^2 - 2\rho x_1 x_2}{2(1-\rho^2)}\right)$$

Another expression of the bivariate Normal copula density is:

$$C(u_1, u_2; \rho) = \int_{0}^{u_1} \Phi \left( \frac{\Phi^{-1}(u_2) - \rho \Phi^{-1}(u_1)}{\sqrt{1 - \rho^2}} \right) du$$

**Student's t copula** In a similar way, the Student's t copula is the dependence function associated with the multivariate Student's t probability distribution:

$$C(u_1, \dots, u_n; \rho, \nu) = T_n \left( T_{\nu}^{-1}(u_1), \dots, T_{\nu}^{-1}(u_n); \rho, \nu \right)$$

By using the definition of the cumulative distribution function:

$$T_n(x_1, \dots, x_n; \rho, \nu) = \int_{-\infty}^{x_1} \dots \int_{-\infty}^{x_n} \frac{\Gamma\left(\frac{\nu+n}{2}\right) |\rho|^{-\frac{1}{2}}}{\Gamma\left(\frac{\nu}{2}\right) (\nu\pi)^{\frac{n}{2}}} \left(1 + \frac{1}{\nu} x^{\top} \rho^{-1} x\right)^{-\frac{\nu+n}{2}} dx$$

we can show that the copula density is then:

$$c(u_1, \dots, u_n; \rho, \nu) = |\rho|^{-\frac{1}{2}} \frac{\Gamma\left(\frac{\nu+n}{2}\right) \left[\Gamma\left(\frac{\nu}{2}\right)\right]^n}{\left[\Gamma\left(\frac{\nu+1}{2}\right)\right]^n \Gamma\left(\frac{\nu}{2}\right)} \frac{\left(1 + \frac{1}{\nu}x^{\top}\rho^{-1}x\right)^{-\frac{\nu+n}{2}}}{\prod_{i=1}^{n} \left(1 + \frac{x_i^2}{\nu}\right)^{-\frac{\nu+1}{2}}}$$

where $x_i = T_{\nu}^{-1}(u_i)$. In the bivariate case, we deduce that the t copula has the following expression:

$$C(u_1, u_2; \rho, \nu) = \int_{-\infty}^{T_{\nu}^{-1}(u_1)} \int_{-\infty}^{T_{\nu}^{-1}(u_2)} \frac{1}{2\pi\sqrt{1-\rho^2}} \left(1 + \frac{x_1^2 + x_2^2 - 2\rho x_1 x_2}{\nu(1-\rho^2)}\right)^{-\frac{\nu+2}{2}} dx_1 dx_2$$

Like the Normal copula, we can obtain another expression, which is easier to manipulate. Let $(X_1, X_2)$ be a random vector whose probability distribution is $T_2(x_1, x_2; \rho, \nu)$. Conditionally to $X_1 = x_1$, we have:

$$\left(\frac{\nu+1}{\nu+x_1^2}\right)^{1/2} \frac{X_2 - \rho x_1}{\sqrt{1-\rho^2}} \sim T_{\nu+1}$$

<footnote n="15">¹⁵In the bivariate case, the parameter $\rho$ is the cross-correlation between $X_1$ and $X_2$, that is the element (1, 2) of the correlation matrix.</footnote>
The conditional distribution $C_{2|1}(u_1, u_2)$ is then equal to:

$$C_{2|1}(u_1, u_2; \rho, \nu) = T_{\nu+1} \left( \left( \frac{\nu+1}{\nu + [T_{\nu}^{-1}(u_1)]^2} \right)^{1/2} \frac{T_{\nu}^{-1}(u_2) - \rho T_{\nu}^{-1}(u_1)}{\sqrt{1 - \rho^2}} \right)$$

We deduce that:

$$C(u_1, u_2; \rho, \nu) = \int_{0}^{u_1} C_{2|1}(u, u_2; \rho, \nu) \, du$$

### A.2.3 Estimation methods

#### Linear regression

Let $Y$ and $X$ be two random vectors. We consider the conditional expectation problem:

$$y = \mathbb{E}[Y | X = x] = m(x)$$

The underlying idea is to find an estimate $\hat{m}(x)$ of the function $m(x)$. In the general case, this problem is extremely difficult to solve. However, if $(Y, X)$ is a Gaussian random vector, the function $m(x)$ can then be determined by considering the Gaussian linear model:

$$Y = \beta^{\top} X + u$$

where $u \sim N(0, \sigma^2)$. Most of the time, the joint distribution of $(Y, X)$ is unknown. In this case, the linear model is estimated by applying least squares techniques<fnref n="16" /> to a given sample $(\mathbf{Y}, \mathbf{X})$:

$$\mathbf{Y} = \mathbf{X}\beta + \mathbf{U}$$

**Derivation of the OLS estimator** We consider a training set of $n$ iid samples $(y_i, x_i)$. For the $i^{th}$ observation, we have:

$$y_i = \sum_{k=1}^{K} \beta_k x_{i,k} + u_i \qquad (A.14)$$

The least squares estimate of the parameter vector $\beta$ is defined as follows:

$$\hat{\beta} = \arg \min \sum_{i=1}^{n} u_{i}^{2}$$

We introduce the following matrix notations: $\mathbf{Y}$ is the $n \times 1$ vector with elements $\mathbf{Y}_i = y_i$, $\mathbf{X}$ is the $n \times K$ matrix defined as follows:

$$\mathbf{X} = \begin{pmatrix} x_{1,1} & & & x_{1,K} \\ & \ddots & & \\ & & \ddots & \\ x_{n,1} & & & x_{n,K} \end{pmatrix}$$

and $\mathbf{U}$ is the $n \times 1$ vector with elements $\mathbf{U}_i = u_i$. In this case, the system of equations (A.14) becomes:

$$\mathbf{Y} = \mathbf{X}\beta + \mathbf{U}$$

<footnote n="16">¹⁶In order to distinguish random variables and observations, we write matrices and vectors that are related to observations in bold style.</footnote>
Let RSS (β) be the residual sum of squares. We have:

$$
\begin{align*}
\text{RSS}(\beta) &= \sum_{i=1}^{n} u_i^2 \\
&= \mathbf{U}^{\top}\mathbf{U} \\
&= \mathbf{Y}^{\top}\mathbf{Y} - 2\beta^{\top}\mathbf{X}^{\top}\mathbf{Y} + \beta^{\top}\mathbf{X}^{\top}\mathbf{X}\beta
\end{align*}
$$

The least squares estimator verifies the set of normal equations $\partial_{\beta}\mathbf{U}^{\top}\mathbf{U} = \mathbf{0}$ and we deduce that $-2\mathbf{X}^{\top}\mathbf{Y} + 2\mathbf{X}^{\top}\mathbf{X}\hat{\beta} = \mathbf{0}$. The expression of the least squares estimator is then:

$$
\hat{\beta} = (\mathbf{X}^{\top}\mathbf{X})^{-1} \mathbf{X}^{\top}\mathbf{Y} \qquad (A.15)
$$

To obtain the expression of $\hat{\beta}$, we only need the assumption that the rank of the matrix $\mathbf{X}$ is $K$. In this case, $\hat{\beta}$ is the solution of the least squares problem. To go further, we assume that $(Y, X)$ is a Gaussian random vector. The solution of the conditional expectation problem $\mathbb{E}[Y | X = x] = m(x)$ is then:

$$
\hat{m}(x) = x^{\top} \hat{\beta} = x^{\top} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{Y}
$$

It means that the prediction of $Y$ given that $X = x$ is equal to $\hat{y} = x^{\top}\hat{\beta}$. If we consider the training data $\mathbf{X}$, we obtain:

$$
\hat{Y} = \hat{m}(\mathbf{X}) = \mathbf{X} (\mathbf{X}^{\top}\mathbf{X})^{-1} \mathbf{X}^{\top}\mathbf{Y} = \mathbf{H}\mathbf{Y}
$$

where $\mathbf{H} = \mathbf{X} (\mathbf{X}^{\top}\mathbf{X})^{-1} \mathbf{X}^{\top}$ is called the ‘hat’ matrix<sup>17</sup>. We notice that $\hat{m}(\mathbf{X})$ is a linear predictor of $\mathbf{Y}$.

**Statistical inference** Because $(Y, X)$ is a Gaussian random vector, it implies that $u = Y - \beta^\top X$ is a Gaussian random variable. We notice that:

$$
\hat{\beta} = (\mathbf{X}^{\top}\mathbf{X})^{-1} \mathbf{X}^{\top}\mathbf{Y} = \beta + (\mathbf{X}^{\top}\mathbf{X})^{-1} \mathbf{X}^{\top}\mathbf{U}
$$

By assuming the exogeneity of the variables $X$ — meaning that $\mathbb{E}[u | X = x] = 0$ — we deduce that $\hat{\beta}$ is an unbiased estimator:

$$
\mathbb{E}[\hat{\beta}] = \beta + (\mathbf{X}^{\top}\mathbf{X})^{-1} \mathbb{E}[\mathbf{X}^{\top}\mathbf{U}] = \beta
$$

$$
\mathbf{H}\mathbf{H} = \mathbf{X} (\mathbf{X}^{\top}\mathbf{X})^{-1} \mathbf{X}^{\top}\mathbf{X} (\mathbf{X}^{\top}\mathbf{X})^{-1} \mathbf{X}^{\top} = \mathbf{X} (\mathbf{X}^{\top}\mathbf{X})^{-1} \mathbf{X}^{\top} = \mathbf{H}
$$

<footnote n="17"><sup>17</sup>We interpret **H** as the orthogonal projection matrix generated by **X** implying that **H** is idempotent, that is **HH** = **H**. Indeed, we have:</footnote>
We recall that $\mathbf{U} \sim \mathcal{N}(0, \sigma^2 I_n)$. It follows that:

$$
\begin{align*}
\operatorname{var}(\hat{\beta}) &= \mathbb{E}\left[(\hat{\beta} - \beta)(\hat{\beta} - \beta)^{\top}\right] \\
&= \mathbb{E}\left[(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{U}\mathbf{U}^{\top}\mathbf{X}(\mathbf{X}^{\top}\mathbf{X})^{-1}\right] \\
&= (\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbb{E}[\mathbf{U}\mathbf{U}^{\top}](\mathbf{X}(\mathbf{X}^{\top}\mathbf{X})^{-1}) \\
&= (\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}(\sigma^2 I_n)\mathbf{X}(\mathbf{X}^{\top}\mathbf{X})^{-1} \\
&= \sigma^2 (\mathbf{X}^{\top}\mathbf{X})^{-1}
\end{align*}
$$

We conclude that:

$$
\hat{\beta} \sim \mathcal{N} \left( \beta, \sigma^2 \left( \mathbf{X}^{\top} \mathbf{X} \right)^{-1} \right)
$$

In most cases, $\sigma^2$ is unknown and we have to estimate it. The vector of residuals is:

$$
\hat{\mathbf{U}} = \mathbf{Y} - \hat{\mathbf{Y}} = \mathbf{Y} - \mathbf{X}\hat{\beta}
$$

We notice that $\mathbb{E}[\hat{\mathbf{U}}] = \mathbf{0}$ and $\operatorname{var}(\hat{\mathbf{U}}) = \sigma^2 (I_n - \mathbf{H})$. Because RSS$(\hat{\beta}) = \hat{\mathbf{U}}^{\top}(I_n - \mathbf{H})\hat{\mathbf{U}}$ is a quadratic form, we can show that:

$$
\hat{\sigma}^2 = \frac{\text{RSS}(\hat{\beta})}{n - K}
$$

is an unbiased estimator of $\sigma^2$ and $\hat{\sigma}^2/\sigma^2 \sim \chi^2_{n-K}$. In order to measure the model quality, we consider the coefficient of determination or $R_c^2$. It is defined as follows:

$$
R_c^2 = 1 - \frac{\text{RSS}(\hat{\beta})}{\text{TSS}}
$$

where $\text{TSS} = \sum_{i=1}^n (y_i - \bar{y})$ is the total sum of squares. We have $R_c^2 \le 1$. A high (resp. low) level indicates a good (resp. bad) goodness-of-fit of the regression model.

Lasso regression

The lasso method consists in adding a $\mathcal{L}_1$-norm penalty function to the optimization function in order to obtain a sparse parameter vector $\theta$:

$$
L_1(\theta) = \|\theta\|_1 = \sum_{k=1}^{K} |\theta_k|
$$

For example, the lasso regression model is specified as follows (Tibshirani, 1996):

$$
y_i = \sum_{k=1}^{K} \beta_k x_{i,k} + u_i \quad \text{s.t.} \quad \sum_{k=1}^{K} |\beta_k| \le \tau
$$

where $\tau$ is a scalar to control the sparsity. Using the notations introduced on page 1169, we have:

$$
\begin{align}
\hat{\beta}(\tau) &= \arg \min (\mathbf{Y} - \mathbf{X}\beta)^{\top} (\mathbf{Y} - \mathbf{X}\beta) \tag{A.16} \\
\text{s.t.} \quad &\|\beta\|_1 \le \tau \nonumber
\end{align}
$$
This problem is equivalent to the Lagrange optimization program $\hat{\beta}(\lambda) = \arg\min \mathcal{L}(\beta; \lambda)$ where<fnref n="18" />:

$$ 
\begin{aligned}
\mathcal{L}(\beta; \lambda) &= \frac{1}{2} (\mathbf{Y} - \mathbf{X}\beta)^{\top} (\mathbf{Y} - \mathbf{X}\beta) + \lambda \|\beta\|_{1} \\
&\propto \frac{1}{2} \beta^{\top} (\mathbf{X}^{\top}\mathbf{X}) \beta - \beta^{\top} (\mathbf{X}^{\top}\mathbf{Y}) + \lambda \|\beta\|_{1}
\end{aligned}
$$

The solution $\hat{\beta}(\lambda)$ can be found by solving the augmented QP program where $\beta = \beta^+ - \beta^-$ under the constraints $\beta^+ \ge 0$ and $\beta^- \ge 0$. We deduce that:

$$ 
\begin{aligned}
\|\beta\|_{1} &= \sum_{k=1}^{K} |\beta_{k}^{+} - \beta_{k}^{-}| \\
&= \sum_{k=1}^{K} |\beta_{k}^{+}| + \sum_{k=1}^{K} |\beta_{k}^{-}| \\
&= \mathbf{1}^{\top}\beta^{+} + \mathbf{1}^{\top}\beta^{-}
\end{aligned}
$$

Since we have:

$$ \beta = ( I_K \quad -I_K ) \begin{pmatrix} \beta^{+} \\ \beta^{-} \end{pmatrix} $$

the augmented QP program is specified as follows:

$$ 
\begin{aligned}
\hat{\theta} &= \arg\min \frac{1}{2} \theta^{\top} Q \theta - \theta^{\top} R \\
\text{s.t.} \quad & \theta \ge 0
\end{aligned}
$$

where $\theta = (\beta^+, \beta^-)$, $\tilde{\mathbf{X}} = (\mathbf{X} - \mathbf{X})$, $Q = \tilde{\mathbf{X}}^\top \tilde{\mathbf{X}}$ and $R = \tilde{\mathbf{X}}^\top \mathbf{Y} - \lambda \cdot \mathbf{1}$. If we denote $A = (I_K - I_K)$, we obtain $\hat{\beta}(\lambda) = A\hat{\theta}$.

**Remark 153** If we consider Problem (A.16), we can also solve it using another augmented QP program:

$$ 
\begin{aligned}
\hat{\theta} &= \arg\min \frac{1}{2} \theta^{\top} Q \theta - \theta^{\top} R \\
\text{s.t.} \quad & \begin{cases} C\theta \ge D \\ \theta \ge 0 \end{cases}
\end{aligned}
$$

where $Q = \tilde{\mathbf{X}}^{\top}\tilde{\mathbf{X}}$, $R = \tilde{\mathbf{X}}^{\top}\mathbf{Y}$, $C = -\mathbf{1}^{\top}$ and $D = -\tau$. We again have $\hat{\beta}(\tau) = A\hat{\theta}$.

We have:

$$ 
\begin{aligned}
\text{RSS}(\beta) &= (\mathbf{Y} - \mathbf{X}\beta)^{\top} (\mathbf{Y} - \mathbf{X}\beta) \\
&= (\mathbf{Y} - \mathbf{X}(\hat{\beta}^{\text{ols}} + \beta - \hat{\beta}^{\text{ols}}))^{\top} (\mathbf{Y} - \mathbf{X}(\hat{\beta}^{\text{ols}} + \beta - \hat{\beta}^{\text{ols}})) \\
&= (\mathbf{Y} - \mathbf{X}\hat{\beta}^{\text{ols}})^{\top} (\mathbf{Y} - \mathbf{X}\hat{\beta}^{\text{ols}}) + 2 (\mathbf{Y} - \mathbf{X}\hat{\beta}^{\text{ols}})^{\top} \mathbf{X} (\beta - \hat{\beta}^{\text{ols}}) + \\
&\quad (\beta - \hat{\beta}^{\text{ols}})^{\top} \mathbf{X}^{\top} \mathbf{X} (\beta - \hat{\beta}^{\text{ols}})
\end{aligned}
$$

<footnote n="18">¹⁸$\tau$ and $\lambda$ are related by the relationship $\tau = \|\hat{\beta}(\lambda)\|_1$.</footnote>
We notice that:

$$ 
\begin{align*}
(*) &= (\mathbf{Y} - \mathbf{X}\hat{\beta}^{\text{ols}})^{\top} \mathbf{X} (\beta - \hat{\beta}^{\text{ols}}) \\
&= (\mathbf{Y}^{\top} - (\hat{\beta}^{\text{ols}})^{\top} \mathbf{X}^{\top}) \mathbf{X} (\beta - \hat{\beta}^{\text{ols}}) \\
&= (\mathbf{Y}^{\top} - ((\mathbf{X}^{\top}\mathbf{X})^{-1} \mathbf{X}^{\top}\mathbf{Y})^{\top} \mathbf{X}^{\top}) \mathbf{X} (\beta - \hat{\beta}^{\text{ols}}) \\
&= (\mathbf{Y}^{\top}\mathbf{X} - ((\mathbf{X}^{\top}\mathbf{X})^{-1} \mathbf{X}^{\top}\mathbf{Y})^{\top} \mathbf{X}^{\top}\mathbf{X}) (\beta - \hat{\beta}^{\text{ols}}) \\
&= (\mathbf{Y}^{\top}\mathbf{X} - \mathbf{Y}^{\top}\mathbf{X}) (\beta - \hat{\beta}^{\text{ols}}) \\
&= 0
\end{align*}
$$

Finally, we obtain:

$$ \text{RSS}(\beta) = \text{RSS}(\hat{\beta}^{\text{ols}}) + (\beta - \hat{\beta}^{\text{ols}})^{\top} \mathbf{X}^{\top}\mathbf{X} (\beta - \hat{\beta}^{\text{ols}}) $$

If we consider the equation RSS ($\beta$) = $c$, we distinguish three cases:

1. if $c < \text{RSS}(\hat{\beta}^{\text{ols}})$, there is no solution;

2. if $c = \text{RSS}(\hat{\beta}^{\text{ols}})$, there is one solution $\beta^* = \hat{\beta}^{\text{ols}}$;

3. if $c > \text{RSS}(\hat{\beta}^{\text{ols}})$, we have:

$$ (\beta - \hat{\beta}^{\text{ols}})^{\top} A (\beta - \hat{\beta}^{\text{ols}}) = 1 $$

where:

$$ A = \frac{\mathbf{X}^{\top}\mathbf{X}}{c - \text{RSS}(\hat{\beta}^{\text{ols}})} $$

The solution $\beta^*$ is an ellipsoid, whose center is $\hat{\beta}^{\text{ols}}$ and principal axes are the eigenvectors of the matrix $A$.

If we add the lasso constraint $\sum_{k=1}^{K} |\beta_k| \le \tau$, the lasso estimator $\hat{\beta}(\tau)$ corresponds to the tangency between the diamond shaped region and the ellipsoid that corresponds to the possible maximum value of $c$. The diamond shape region due to the lasso constraint ensures that the lasso estimator is sparse:

$$ \exists \eta > 0 : \forall \tau < \eta, \min(\hat{\beta}_1(\tau), \dots, \hat{\beta}_K(\tau)) = 0 $$

For example, the two-dimensional case is represented in Figure A.7. We notice that $\hat{\beta}_1(\tau)$ is equal to zero if $\tau < \eta$. This sparsity property is central for understanding the variable selection procedure.

**Example 69** Using the data given in Table A.2, we consider the linear regression model:

$$ y_i = \beta_0' + \sum_{k=1}^{5} \beta_k' x_{i,k} + u_i \quad (A.17) $$

The objective is to determine the importance of each variable.
<figure><img src="image_539.png" /><figcaption>Figure A.7: Interpretation of the lasso regression</figcaption></figure>

The lasso method can be used for ranking the variables. For that, we consider the following linear regression:

$$ \tilde{y}_i = \sum_{k=1}^{5} \beta_k \tilde{x}_{i,k} + u_i $$

where $\tilde{y}_i$ and $\tilde{x}_{i,k}$ are the standardized data<fnref n="19" />:

$$ \frac{y_i - \bar{y}}{s_y} = \sum_{k=1}^{5} \beta_k \left( \frac{x_{i,k} - \bar{x}_k}{s_{x_k}} \right) + u_i \quad (A.18) $$

Linear regressions (A.17) and (A.18) are related by the following equation:

$$ y_i = \left( \bar{y} - \sum_{k=1}^{5} \frac{s_y \beta_k}{s_{x_k}} \bar{x}_k \right) + \sum_{k=1}^{5} \frac{s_y \beta_k}{s_{x_k}} x_{i,k} + s_y u_i $$

We deduce that $\beta'_0 = \bar{y} - \sum_{k=1}^{5} (s_y/s_{x_k}) \beta_k \bar{x}_k$ and $\beta'_k = (s_y/s_{x_k}) \beta_k$. When performing lasso regression, we always standardize the data in order to obtain comparable beta's. Otherwise, the penalty function $||\beta||_1$ does not make a lot of sense. In Table A.3, we have estimated the lasso coefficients $\beta_k(\lambda)$ for different values of the shrinkage parameter $\lambda$. When $\lambda = 0$, we obtain the OLS estimate, and the lasso regression selects all the available variables. When $\lambda \to \infty$, the solution is $\hat{\beta}(\infty) = 0$, and the lasso regression selects no explanatory variables. In Table A.3, we verify that the number of selected variables is a decreasing function of $\lambda$. For instance, the lasso regression selects respectively four and three variables when $\lambda$ is equal to 0.9 and 2.5. It follows that the most important variable is the third one, followed by the first, second, fourth and fifth variables.

<footnote n="19">¹⁹The notations $\bar{x}_k$ and $s_{x_k}$ represent the mean and the standard deviation of the data $\{x_{i,k}, i = 1, \dots, n\}$.</footnote>
<table><caption>Table A.2: Data of the lasso regression problem</caption><thead><tr><th><i>i</i></th><th><i>y</i></th><th><i>x</i><sub>1</sub></th><th><i>x</i><sub>2</sub></th><th><i>x</i><sub>3</sub></th><th><i>x</i><sub>4</sub></th><th><i>x</i><sub>5</sub></th></tr></thead><tbody><tr><td>1</td><td>3.1</td><td>2.8</td><td>4.3</td><td>0.3</td><td>2.2</td><td>3.5</td></tr><tr><td>2</td><td>24.9</td><td>5.9</td><td>3.6</td><td>3.2</td><td>0.7</td><td>6.4</td></tr><tr><td>3</td><td>27.3</td><td>6.0</td><td>9.6</td><td>7.6</td><td>9.5</td><td>0.9</td></tr><tr><td>4</td><td>25.4</td><td>8.4</td><td>5.4</td><td>1.8</td><td>1.0</td><td>7.1</td></tr><tr><td>5</td><td>46.1</td><td>5.2</td><td>7.6</td><td>8.3</td><td>0.6</td><td>4.5</td></tr><tr><td>6</td><td>45.7</td><td>6.0</td><td>7.0</td><td>9.6</td><td>0.6</td><td>0.6</td></tr><tr><td>7</td><td>47.4</td><td>6.1</td><td>1.0</td><td>8.5</td><td>9.6</td><td>8.6</td></tr><tr><td>8</td><td>-1.8</td><td>1.2</td><td>9.6</td><td>2.7</td><td>4.8</td><td>5.8</td></tr><tr><td>9</td><td>20.8</td><td>3.2</td><td>5.0</td><td>4.2</td><td>2.7</td><td>3.6</td></tr><tr><td>10</td><td>6.8</td><td>0.5</td><td>9.2</td><td>6.9</td><td>9.3</td><td>0.7</td></tr><tr><td>11</td><td>12.9</td><td>7.9</td><td>9.1</td><td>1.0</td><td>5.9</td><td>5.4</td></tr><tr><td>12</td><td>37.0</td><td>1.8</td><td>1.3</td><td>9.2</td><td>6.1</td><td>8.3</td></tr><tr><td>13</td><td>14.7</td><td>7.4</td><td>5.6</td><td>0.9</td><td>5.6</td><td>3.9</td></tr><tr><td>14</td><td>-3.2</td><td>2.3</td><td>6.6</td><td>0.0</td><td>3.6</td><td>6.4</td></tr><tr><td>15</td><td>44.3</td><td>7.7</td><td>2.2</td><td>6.5</td><td>1.3</td><td>0.7</td></tr></tbody></table>

<table><caption>Table A.3: Results of the lasso regression</caption><thead><tr><th>λ</th><th>0.0</th><th>0.9</th><th>2.5</th><th>5.5</th><th>7.5</th></tr></thead><tbody><tr><td>β̂<sub>1</sub>(λ)</td><td>0.4586</td><td>0.4022</td><td>0.3163</td><td>0.1130</td><td></td></tr><tr><td>β̂<sub>2</sub>(λ)</td><td>-0.1849</td><td>-0.2005</td><td>-0.1411</td><td></td><td></td></tr><tr><td>β̂<sub>3</sub>(λ)</td><td>0.8336</td><td>0.7265</td><td>0.5953</td><td>0.3951</td><td>0.2462</td></tr><tr><td>β̂<sub>4</sub>(λ)</td><td>-0.1893</td><td>-0.1102</td><td></td><td></td><td></td></tr><tr><td>β̂<sub>5</sub>(λ)</td><td>0.0931</td><td></td><td></td><td></td><td></td></tr><tr><td>||β̂(λ)||<sub>1</sub></td><td>1.7595</td><td>1.4395</td><td>1.0527</td><td>0.5081</td><td>0.2462</td></tr><tr><td>RSS(β̂(λ))</td><td>0.0118</td><td>0.0304</td><td>0.1180</td><td>0.4076</td><td>0.6306</td></tr><tr><td>R<sup>2</sup><sub>c</sub></td><td>0.9874</td><td>0.9674</td><td>0.8735</td><td>0.5633</td><td>0.3244</td></tr><tr><td>df(model)</td><td>5</td><td>4</td><td>3</td><td>2</td><td>1</td></tr></tbody></table>

In Figure A.8, we have reported the path of the lasso estimate $\hat{\beta}(\lambda)$ with respect to the scaling factor $\tau^* \in [0, 1]$, which is defined as follows:

$$ \tau^* = \frac{\tau}{\tau_{\max}} = \frac{\|\hat{\beta}(\lambda)\|_1}{\|\hat{\beta}(0)\|_1} $$

$\tau^*$ is equal to zero when $\lambda \to \infty$ (no selected variable) and one when $\lambda = 0$, which corresponds to the OLS case. From this path, we verify the lasso ordering $x_3 \succ x_1 \succ x_2 \succ x_4 \succ x_5$.

## State space models

A state space model (SSM) includes a measurement equation and a transition equation. In the measurement equation, we define the relationship between an observable system and state variables, whereas the transition equation describes the dynamics of state variables. Generally, the state vector
<figure><img src="image_540.png" /><figcaption>Figure A.8: Variable selection with the lasso regression</figcaption></figure>

$\alpha_t$ is generated by a Markov linear process<fnref n="20" />:

$$\alpha_t = T_t \alpha_{t-1} + c_t + R_t \eta_t$$

where $\alpha_t$ is a $m \times 1$ vector, $T_t$ is a $m \times m$ matrix, $c_t$ is a $m \times 1$ vector and $R_t$ is a $m \times p$ matrix. In the case of a linear SSM, the measurement equation is given by:

$$y_t = Z_t \alpha_t + d_t + \epsilon_t$$

where $y_t$ is a $n$-dimensional time series, $Z_t$ is a $n \times m$ matrix, $d_t$ is a $n \times 1$ vector. We also assume that $\eta_t$ and $\epsilon_t$ are two independent white noise processes of dimension $p$ and $n$ with covariance matrices $Q_t$ and $H_t$.

**Kalman filtering** In the state space model, the variable $y_t$ is observable, but it is generally not the case of the state vector $\alpha_t$. The Kalman filter is a statistical tool to estimate the distribution function of $\alpha_t$. Let $\alpha_0 \sim N(\hat{\alpha}_0, P_0)$ the initial position of the state vector. We note $\hat{\alpha}_{t|t}$ (or $\hat{\alpha}_t$) and $\hat{\alpha}_{t|t-1}$ the optimal estimators of $\alpha_t$ given the available information until time $t$ and $t-1$:

$$\hat{\alpha}_{t|t} = \mathbb{E}[\alpha_t | \mathcal{F}_t] \\ \hat{\alpha}_{t|t-1} = \mathbb{E}[\alpha_t | \mathcal{F}_{t-1}]$$

$P_{t|t}$ (or $P_t$) and $P_{t|t-1}$ are the covariance matrices associated to $\hat{\alpha}_{t|t}$ and $\hat{\alpha}_{t|t-1}$:

$$P_{t|t} = \mathbb{E}[(\hat{\alpha}_{t|t} - \alpha_t)(\hat{\alpha}_{t|t} - \alpha_t)^{\top}] \\ P_{t|t-1} = \mathbb{E}[(\hat{\alpha}_{t|t-1} - \alpha_t)(\hat{\alpha}_{t|t-1} - \alpha_t)^{\top}]$$

<footnote n="20">²⁰The presentation is based on the book of Harvey (1990).</footnote>
These different quantities are calculated thanks to the Kalman filter, which consists in a recursive algorithm<sup>21</sup> (Harvey, 1990):

$$
\left\{
\begin{aligned}
\hat{\alpha}_{t|t-1} &= T_t \hat{\alpha}_{t-1|t-1} + c_t \\
P_{t|t-1} &= T_t P_{t-1|t-1} T_t^\top + R_t Q_t R_t^\top \\
\hat{y}_{t|t-1} &= Z_t \hat{\alpha}_{t|t-1} + d_t \\
v_t &= y_t - \hat{y}_{t|t-1} \\
F_t &= Z_t P_{t|t-1} Z_t^\top + H_t \\
\hat{\alpha}_{t|t} &= \hat{\alpha}_{t|t-1} + P_{t|t-1} Z_t^\top F_t^{-1} v_t \\
P_{t|t} &= (I_m - P_{t|t-1} Z_t^\top F_t^{-1} Z_t) P_{t|t-1}
\end{aligned}
\right.
$$

where $\hat{y}_{t|t-1} = \mathbb{E}[y_t | \mathcal{F}_{t-1}]$ is the best estimator of $y_t$ given the available information until time $t-1$, $v_t$ is the innovation process and $F_t$ is the associated covariance matrix.

**Remark 154** Harvey (1990) showed that we can directly calculate $\hat{\alpha}_{t+1|t}$ from $\hat{\alpha}_{t|t-1}$:

$$ \hat{\alpha}_{t+1|t} = (T_{t+1} - K_t Z_t) \hat{\alpha}_{t|t-1} + K_t y_t + c_{t+1} - K_t d_t $$

where $K_t = T_{t+1}P_{t|t-1}Z_t^\top F_t^{-1}$ is the gain matrix. It follows that:

$$ \hat{\alpha}_{t+1|t} = T_{t+1}\hat{\alpha}_{t|t-1} + c_{t+1} + K_t (y_t - Z_t\hat{\alpha}_{t|t-1} - d_t) $$

By recognizing the innovation process $v_t$, we obtain the following innovation representation:

$$
\left\{
\begin{aligned}
y_t &= Z_t \hat{\alpha}_{t|t-1} + d_t + v_t \\
\hat{\alpha}_{t+1|t} &= T_{t+1} \hat{\alpha}_{t|t-1} + c_{t+1} + K_t v_t
\end{aligned}
\right.
$$

**Estimation of unknown parameters** In many cases, the state space model depends on certain parameters that are unknown. Given a set $\theta$ of values for these unknown parameters, the Kalman filter may be applied to estimate the state vector $\alpha_t$. We have:

$$ v_t \sim \mathcal{N}(0, F_t) $$

where $v_t = y_t - \hat{y}_{t|t-1}$ is the innovation at time $t$ and $F_t = Z_t P_{t|t-1} Z_t^\top + H_t$ is the covariance matrix. If we change $\theta$ and we run the Kalman filter, we will obtain other values of $v_t$ and $F_t$, meaning that $v_t$ and $F_t$ depend on $\theta$. This is why we can write $v_t(\theta)$ and $F_t(\theta)$. We deduce that the likelihood function of the sample $\{y_1, \dots, y_T\}$ is equal to:

$$ \ell(\theta) = -\frac{nT}{2} \ln(2\pi) - \frac{1}{2} \sum_{t=1}^{T} \left( \ln|F_t(\theta)| + v_t(\theta)^{\top} F_t(\theta)^{-1} v_t(\theta) \right) $$

We can then estimate the vector $\theta$ of unknown parameters by the method of maximum likelihood:

$$ \hat{\theta} = \arg\max \ell(\theta) $$

Once the ML estimate $\hat{\theta}$ is found, we can run again<sup>22</sup> the Kalman filter to estimate the other quantities $\hat{\alpha}_{t|t-1}$, $\hat{\alpha}_{t|t}$, $P_{t|t-1}$ and $P_{t|t}$.

<footnote n="21"><sup>21</sup>The algorithm is initialized with values $\hat{\alpha}_{0|0} = \hat{\alpha}_0$ and $P_{0|0} = P_0$.</footnote>

<footnote n="22"><sup>22</sup>We say again, because computing the log-likelihood function requires one Kalman filter run, implying that many Kalman filter runs are used for maximizing the log-likelihood function.</footnote>
**Time-invariant state space model** We consider the time-invariant model:

$$ egin{cases} y_t = Z\alpha_t + d + \epsilon_t \\ \alpha_t = T\alpha_{t-1} + c + R\eta_t \end{cases} $$

where $\epsilon_t \sim \mathcal{N}(0, H)$ and $\eta_t \sim \mathcal{N}(0, Q)$. If the state space model converges to a steady state, the estimators $(\hat{\alpha}_{\infty}, P_{\infty})$ must satisfy the following equations:

$$ egin{cases} \hat{\alpha}_{\infty} = T\hat{\alpha}_{\infty} + c \\ P_{\infty} = T P_{\infty} T^{\top} + RQR^{\top} \end{cases} $$

It follows that the solution is:

$$ egin{cases} \hat{\alpha}_{\infty} = (I_m - T)^{-1} c \\ \mathrm{vec}(P_{\infty}) = (I_{m^2} - T \otimes T)^{-1} \mathrm{vec}(RQR^{\top}) \end{cases} $$

where $\hat{\alpha}_{\infty}$ and $P_{\infty}$ are the unconditional mean and covariance matrix of $\alpha_t$. Without any knowledge of the initial position $\alpha_0$, the best way to define $\hat{\alpha}_0$ and $P_0$ is then to use the steady state:

$$ egin{cases} \hat{\alpha}_0 = \hat{\alpha}_{\infty} \\ P_0 = P_{\infty} \end{cases} $$

In many state space models, the matrices $T$, $c$, $R$ and $Q$ depend on unknown parameters $\theta$, implying that $\hat{\alpha}_{\infty}$ and $P_{\infty}$ also depend on $\theta$. This means that when maximizing the log-likelihood function, the Kalman filter is initialized by values of $\hat{\alpha}_0$ and $P_0$ that depend on $\theta$. This is the main difference with time-varying state space model since the Kalman filter is initialized by fixed values of $\hat{\alpha}_0$ and $P_0$.

## A.3 Stochastic analysis

### A.3.1 Stochastic optimal control

### A.3.2 Jump-diffusion processes

### A.4 Spatial data

#### A.4.1 Spherical coordinates

#### A.4.2 Geographic coordinate systems

#### A.4.3 Network common data form
# Appendix B

## Solutions to the Tutorial Exercises

### B.1 Exercises related to ESG risk

#### B.1.1 Score normalization when the features are independent (Exercise 2.4.1)

#### B.1.2 Probability distribution of ESG scores (Exercise 2.4.4)

1. (a) We have:

$$ \mathbf{S}_i^{(\text{ESG})} = 0.4 \times \mathbf{S}_i^{(\text{E})} + 0.4 \times \mathbf{S}_i^{(\text{S})} + 0.2 \times \mathbf{S}_i^{(\text{G})} $$

We deduce the following results:

<table><thead><tr><th>Issuer</th><th>#1</th><th>#2</th><th>#3</th><th>#4</th><th>#5</th><th>#6</th><th>#7</th><th>#8</th></tr></thead><tbody><tr><td><b>S</b><sub>i</sub><sup>(E)</sup></td><td>-2.80</td><td>-1.80</td><td>-1.75</td><td>0.60</td><td>0.75</td><td>1.30</td><td>1.90</td><td>2.70</td></tr><tr><td><b>S</b><sub>i</sub><sup>(S)</sup></td><td>-1.70</td><td>-1.90</td><td>0.75</td><td>-1.60</td><td>1.85</td><td>1.05</td><td>0.90</td><td>0.70</td></tr><tr><td><b>S</b><sub>i</sub><sup>(G)</sup></td><td>0.30</td><td>-0.70</td><td>-2.75</td><td>2.60</td><td>0.45</td><td>2.35</td><td>2.20</td><td>1.70</td></tr><tr><td><b>S</b><sub>i</sub><sup>(ESG)</sup></td><td>-1.74</td><td>-1.62</td><td>-0.95</td><td>0.12</td><td>1.13</td><td>1.41</td><td>1.56</td><td>1.70</td></tr></tbody></table>

(b) We obtain:

$$ \mathbf{S}^{(\text{ESG})}(x_{\text{ew}}) = \sum_{i=1}^{8} x_{\text{ew},i} \times \mathbf{S}_{i}^{(\text{ESG})} = 0.2013 $$

2. (a) We have:

$$ \mathbf{S} \left( x_{\text{ew}}^{(n)} \right) = \sum_{i=1}^{n} x_{\text{ew},i}^{(n)} \times \mathbf{S}_i = \frac{1}{n} \sum_{i=1}^{n} \mathbf{S}_i $$

We deduce that $\mathbf{S} \left( x_{\text{ew}}^{(n)} \right)$ follows a Gaussian distribution. Its mean is equal to:

$$ \mathbb{E} \left[ \mathbf{S} \left( x_{\text{ew}}^{(n)} \right) \right] = \frac{1}{n} \sum_{i=1}^{n} \mathbb{E} [\mathbf{S}_i] = 0 $$

whereas its standard deviation is equal to:

$$ \sigma \left( \mathbf{S} \left( x_{\text{ew}}^{(n)} \right) \right) = \sqrt{\frac{1}{n^2} \sum_{i=1}^{n} \sigma^2 (\mathbf{S}_i)} = \frac{1}{\sqrt{n}} $$
Finally, we deduce that:

$$ \mathcal{S}(x_{\text{ew}}^{(n)}) \sim \mathcal{N}\left(0, \frac{1}{n}\right) $$

(b) The behavior of a well-diversified portfolio is close to an equally-weighted portfolio with $n$ sufficiently large. Therefore, the ESG score is close to zero because we have:

$$ \lim_{n \to \infty} \mathcal{S}(x_{\text{ew}}^{(n)}) = 0 $$

(c) We have:

$$ f_{\alpha}(t) = \alpha t^{\alpha-1} $$

The probability density function $f_{\alpha}(t)$ is reported in Figure B.1. We notice that the function $f_{\alpha}(t)$ tends to the dirac delta function when $\alpha$ tends to infinity:

$$ \lim_{\alpha \to \infty} f_{\alpha}(t) = \delta_{1}(t) = \begin{cases} 0 & \text{if } t \neq 1 \\ +\infty & \text{if } t = 1 \end{cases} $$

<figure><img src="image_541.png" /><figcaption>Figure B.1: Probability density function $f_{\alpha}(t)$</figcaption></figure>

(d) To simulate $T_i$, we use the property of the probability integral transform: $U_i = \mathbf{F}_{\alpha}(T_i) \sim \mathcal{U}_{[0,1]}$. We deduce that:

$$ T_i = \mathbf{F}_{\alpha}^{-1}(U_i) = U_i^{1/\alpha} $$

The algorithm for simulating the portfolio $x$ is then the following:

i. We simulate $n$ independent uniform random numbers $(u_1, \dots, u_n)$;
ii. We compute the random variates $(t_1, \dots, t_n)$ where:

$$t_i = u_i^{1/\alpha}$$

iii. We calculate the normalization constant:

$$c = \left( \sum_{i=1}^{n} t_i \right)^{-1} = \left( \sum_{i=1}^{n} u_i^{1/\alpha} \right)^{-1}$$

iv. We deduce the portfolio weights $x = (x_1, \dots, x_n)$:

$$x_i = c \cdot t_i = c \cdot u_i^{1/\alpha} = \frac{u_i^{1/\alpha}}{\sum_{j=1}^{n} u_j^{1/\alpha}}$$

In Figure B.2, we have represented the composition of the portfolio $x$ for the 4 values of $\alpha$. The weights are ranked in descending order. We deduce that the portfolio $x$ is uniform when $\alpha \to \infty$. The parameter $\alpha$ controls the concentration of the portfolio. Indeed, when $\alpha$ is small, the portfolio is highly concentrated. It follows that the Herfindahl index $H_\alpha(x)$ of the portfolio weights is a decreasing function of the parameter $\alpha$.

<figure><img src="image_542.png" /><figcaption>Figure B.2: Repartition of the portfolio weights in descending order</figcaption></figure>

(e) We simulate $x = (x_1, \dots, x_n)$ using the previous algorithm. The vector of ESG scores $\mathbf{S} = (\mathbf{S}_1, \dots, \mathbf{S}_n)$ is generated with normally-distributed random variables since we have $\mathbf{S}_i \sim \mathcal{N}(0, 1)$. We deduce that the simulated value of the portfolio ESG score $\mathbf{S}(x)$ is equal to:

$$\mathbf{S}(x) = \sum_{i=1}^{n} x_i \cdot \mathbf{S}_i$$
We replicate the simulation of $S(x)$ 50 000 times and draw the corresponding histogram in Figure B.3. We also report the fitted Gaussian distribution. We observe that the portfolio's ESG score $S(x)$ is equal to zero on average, and its variance is an increasing function of the portfolio concentration.

<figure><img src="image_543.png" /><figcaption>Figure B.3: Histogram of the portfolio ESG score $S(x)$</figcaption></figure>

(f) Since $x_i \sim cT_i$, $x_i$ is an increasing function of $T_i$. We deduce that the copula function of $(T_i, S_i)$ is the same as the copula function of $(x_i, S_i)$. To simulate the Normal copula function $C(u, v)$, we use the transformation algorithm based on the Cholesky decomposition:

$$ \begin{cases} u_i = \Phi(g'_i) \\ v_i = \Phi(\rho g'_i + \sqrt{1 - \rho^2} g''_i) \end{cases} $$

where $g'_i$ and $g''_i$ are two independent random numbers from the probability distribution $N(0, 1)$. Here is the algorithm to simulate the portfolios' ESG score $S(x)$:

i. We simulate $n$ independent normally-distributed random numbers $g'_i$ and $g''_i$ and compute $(u_i, v_i)$:

$$ \begin{cases} u_i = \Phi(g'_i) \\ v_i = \Phi(\rho g'_i + \sqrt{1 - \rho^2} g''_i) \end{cases} $$

ii. We compute the random variates $(t_1, \dots, t_n)$ where $t_i = u_i^{1/\alpha}$;

iii. We deduce the vector of weights $x = (x_1, \dots, x_n)$:

$$ x_i = t_i / \sum_{j=1}^{n} t_j $$
iv. We simulate the vector of scores $\mathbf{S} = (\mathbf{S}_1, \dots, \mathbf{S}_n)$:

$$ \mathbf{S}_i = \Phi^{-1}(v_i) = \rho g'_i + \sqrt{1 - \rho^2} g''_i $$

v. We calculate the portfolio score:

$$ \mathbf{S}(x) = \sum_{i=1}^{n} x_i \cdot \mathbf{S}_i $$

<figure><img src="image_544.png" /><figcaption>Figure B.4: Histogram of the portfolio ESG score $\mathbf{S}(x)$ ($\rho = 50\%$)</figcaption></figure>

(g) In the independent case, we found that $\mathbb{E}[\mathbf{S}(x)] = 0$. In Figure B.4, we notice that $\mathbb{E}[\mathbf{S}(x)] \neq 0$ when $\rho$ is equal to 50%. Indeed, we obtain:

$$ \mathbb{E}[\mathbf{S}(x)] = \begin{cases} 0.418 & \text{if } \alpha = 0.5 \\ 0.210 & \text{if } \alpha = 1.5 \\ 0.142 & \text{if } \alpha = 2.5 \\ 0.006 & \text{if } \alpha = 70.0 \end{cases} $$

(h) In Figure B.5, we notice that there is a positive relationship between $\rho$ and $\mathbb{E}[\mathbf{S}(x)]$ and the slope increases with the concentration of the portfolio.

(i) Big cap companies have more (financial and human) resources to develop an ESG policy than small cap companies. Therefore, we observe a positive correlation between the market capitalization and the ESG score of an issuer. It follows that ESG portfolios have generally a size bias. For instance, we generally observe that cap-weighted indexes have an ESG score which is greater than the average of ESG scores. In the previous questions, we verify that $\mathbb{E}[\mathbf{S}(x)] \ge \mathbb{E}[\mathbf{S}]$ when the Herfindahl index of the portfolio $x$ is high and the correlation between $x_i$ and $\mathbf{S}_i$ is positive.
<figure><img src="image_545.png" /><figcaption>Figure B.5: Relationship between $\rho$ and $\mathbb{E}[\mathbf{S}(x)]$</figcaption></figure>

3. (a) We have:

$$ 
\begin{aligned}
\mathrm{Pr}\{\mathcal{R} = \mathbf{A}\} &= \mathrm{Pr}\{\mathbf{S} \geq 1.5\} \\
&= 1 - \Phi(1.5) \\
&= 6.68\%
\end{aligned}
$$

and:

$$ 
\begin{aligned}
\mathrm{Pr}\{\mathcal{R} = \mathbf{B}\} &= \mathrm{Pr}\{0 \leq \mathbf{S} < 1.5\} \\
&= \Phi(1.5) - \Phi(0) \\
&= 43.32\%
\end{aligned}
$$

Since the Gaussian distribution is symmetric around 0, we also have:

$$ \mathrm{Pr}\{\mathcal{R} = \mathbf{C}\} = \mathrm{Pr}\{\mathcal{R} = \mathbf{B}\} = 43.32\% $$

and:

$$ \mathrm{Pr}\{\mathcal{R} = \mathbf{D}\} = \mathrm{Pr}\{\mathcal{R} = \mathbf{A}\} = 6.68\% $$

The mapping function is then equal to:

$$ \mathcal{M}_{\text{appring}}(\mathbf{S}) = \begin{cases} \mathbf{A} & \text{if } \mathbf{S} < -1.5 \\ \mathbf{B} & \text{if } -1.5 \leq \mathbf{S} < 0 \\ \mathbf{C} & \text{if } 0 \leq \mathbf{S} < 1.5 \\ \mathbf{D} & \text{if } \mathbf{S} \geq 1.5 \end{cases} $$
(b) Since we have:

$$ \mathrm{Pr}\{\mathcal{R}(t) = \mathbf{A}\} = \mathrm{Pr}\{\mathcal{R}(t) = \mathbf{B}\} = \mathrm{Pr}\{\mathcal{R}(t) = \mathbf{C}\} = \mathrm{Pr}\{\mathcal{R}(t) = \mathbf{D}\} $$

and:

$$ \mathrm{Pr}\{\mathcal{R}(t) = \mathbf{A}\} + \mathrm{Pr}\{\mathcal{R}(t) = \mathbf{B}\} + \mathrm{Pr}\{\mathcal{R}(t) = \mathbf{C}\} + \mathrm{Pr}\{\mathcal{R}(t) = \mathbf{D}\} = 1 $$

we deduce that:

$$ \mathrm{Pr}\{\mathcal{R}(t) = \mathbf{A}\} = \frac{1}{4} = 25\% $$

and $\mathrm{Pr}\{\mathcal{R}(t) = \mathbf{B}\} = \mathrm{Pr}\{\mathcal{R}(t) = \mathbf{C}\} = \mathrm{Pr}\{\mathcal{R}(t) = \mathbf{D}\} = 25\%$. We want to find the breakpoints $(s_1, s_2, s_3)$ such that:

$$ \left\{ \begin{array}{l} \mathrm{Pr}\{\mathbf{S} < s_1\} = 25\% \\ \mathrm{Pr}\{s_1 \le \mathbf{S} < s_2\} = 25\% \\ \mathrm{Pr}\{s_2 \le \mathbf{S} < s_3\} = 25\% \\ \mathrm{Pr}\{\mathbf{S} \ge s_3\} = 25\% \end{array} \right. $$

We deduce that:

$$ \left\{ \begin{array}{l} s_1 = \Phi^{-1}(0.25) = -0.6745 \\ s_2 = \Phi^{-1}(0.50) = 0 \\ s_3 = \Phi^{-1}(0.75) = +0.6745 \end{array} \right. $$

The mapping function is then given by:

$$ \mathcal{M}_{\text{appring}}(\mathbf{S}) = \begin{cases} \mathbf{A} & \text{if } \mathbf{S} < -0.6745 \\ \mathbf{B} & \text{if } -0.6745 \le \mathbf{S} < 0 \\ \mathbf{C} & \text{if } 0 \le \mathbf{S} < 0.6745 \\ \mathbf{D} & \text{if } \mathbf{S} \ge 0.6745 \end{cases} $$

(c) We have:

$$ \left\{ \begin{array}{l} s_1 = \Phi^{-1}(0.10) = -1.2816 \\ s_2 = \Phi^{-1}(0.50) = 0 \\ s_3 = \Phi^{-1}(0.90) = +1.2816 \end{array} \right. $$

We deduce that the mapping function is equal to:

$$ \mathcal{M}_{\text{appring}}(\mathbf{S}) = \begin{cases} \mathbf{A} & \text{if } \mathbf{S} < -1.2816 \\ \mathbf{B} & \text{if } -1.2816 \le \mathbf{S} < 0 \\ \mathbf{C} & \text{if } 0 \le \mathbf{S} < 1.2816 \\ \mathbf{D} & \text{if } \mathbf{S} \ge 1.2816 \end{cases} $$

4. (a) The joint distribution of $(\mathbf{S}(t-1), \Delta\mathbf{S}(t))$ is:

$$ \begin{pmatrix} \mathbf{S}(t-1) \\ \Delta\mathbf{S}(t) \end{pmatrix} \sim \mathcal{N} \left( \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 1 & 0 \\ 0 & \sigma^2 \end{pmatrix} \right) $$

(b) Since we have:

$$ \mathbf{S}(t) = \mathbf{S}(t-1) + \Delta\mathbf{S}(t) $$

we deduce that:

$$ \begin{pmatrix} \mathbf{S}(t-1) \\ \mathbf{S}(t) \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 1 & 1 \end{pmatrix} \begin{pmatrix} \mathbf{S}(t-1) \\ \Delta\mathbf{S}(t) \end{pmatrix} $$
We conclude that $(\mathbf{S}(t-1), \mathbf{S}(t))$ is a Gaussian random vector. We have:

$$ \mathrm{var}(\mathbf{S}(t)) = 1 + \sigma^2 $$

and:

$$ \begin{aligned} \mathrm{cov}(\mathbf{S}(t-1), \mathbf{S}(t)) &= \mathbb{E}[\mathbf{S}(t-1) \cdot \mathbf{S}(t)] \\ &= \mathbb{E}[\mathbf{S}^2(t-1) + \mathbf{S}(t-1) \cdot \Delta\mathbf{S}(t)] \\ &= 1 \end{aligned} $$

It follows that:

$$ \begin{pmatrix} \mathbf{S}(t-1) \\ \mathbf{S}(t) \end{pmatrix} \sim \mathcal{N}(\mathbf{0}_2, \Sigma_\sigma) $$

where $\Sigma_\sigma$ is the covariance matrix:

$$ \Sigma_\sigma = \begin{pmatrix} 1 & 1 \\ 1 & 1 + \sigma^2 \end{pmatrix} $$

(c) We have:

$$ \mathrm{Pr}\{\mathcal{R}(t-1) = \mathcal{R}_k\} = \mathrm{Pr}\{s_{k-1} \le \mathbf{S}(t-1) < s_k\} = \Phi(s_k) - \Phi(s_{k-1}) $$

(d) We have:

$$ \begin{aligned} (*) &= \mathrm{Pr}\{\mathcal{R}(t) = \mathcal{R}_k, \mathcal{R}(t-1) = \mathcal{R}_j\} \\ &= \mathrm{Pr}\{s_{k-1} \le \mathbf{S}(t) < s_k, s_{j-1} \le \mathbf{S}(t-1) < s_j\} \\ &= \Phi_2(s_j, s_k; \Sigma_\sigma) - \Phi_2(s_{j-1}, s_k; \Sigma_\sigma) - \Phi_2(s_j, s_{k-1}; \Sigma_\sigma) + \Phi_2(s_{j-1}, s_{k-1}; \Sigma_\sigma) \end{aligned} $$

where $\Phi_2(x, y; \Sigma_\sigma)$ is the bivariate Normal cdf with covariance matrix $\Sigma_\sigma$.

(e) We have:

$$ \begin{aligned} p_{j,k} &= \mathrm{Pr}\{\mathcal{R}(t) = \mathcal{R}_k | \mathcal{R}(t-1) = \mathcal{R}_j\} \\ &= \frac{\mathrm{Pr}\{\mathcal{R}(t) = \mathcal{R}_k, \mathcal{R}(t-1) = \mathcal{R}_j\}}{\mathrm{Pr}\{\mathcal{R}(t-1) = \mathcal{R}_j\}} \\ &= \frac{\Phi_2(s_j, s_k; \Sigma_\sigma) + \Phi_2(s_{j-1}, s_{k-1}; \Sigma_\sigma)}{\Phi_2(s_j, s_k; \Sigma_\sigma) + \Phi_2(s_{j-1}, s_{k-1}; \Sigma_\sigma)} - \frac{\Phi_2(s_j, s_k; \Sigma_\sigma) + \Phi_2(s_{j-1}, s_{k-1}; \Sigma_\sigma)}{\Phi_2(s_j, s_k; \Sigma_\sigma) + \Phi_2(s_{j-1}, s_{k-1}; \Sigma_\sigma)} \end{aligned} $$

(f) We have:

$$ \begin{aligned} \mathcal{T}(\mathcal{R}_k) &= \mathrm{Pr}\{\mathcal{R}(t) \neq \mathcal{R}_k | \mathcal{R}(t-1) = \mathcal{R}_k\} \\ &= 1 - \mathrm{Pr}\{\mathcal{R}(t) = \mathcal{R}_k | \mathcal{R}(t-1) = \mathcal{R}_k\} \\ &= 1 - p_{k,k} \end{aligned} $$

(g) We have:

$$ \begin{aligned} \mathcal{T}(\mathcal{R}_1, \dots, \mathcal{R}_K) &= \sum_{k=1}^{K} \mathrm{Pr}\{\mathcal{R}(t-1) = \mathcal{R}_k\} \cdot \mathcal{T}(\mathcal{R}_k) \\ &= \sum_{k=1}^{K} \mathrm{Pr}\{\mathcal{R}(t) \neq \mathcal{R}_k, \mathcal{R}(t-1) = \mathcal{R}_k\} \end{aligned} $$
<table><caption>Table B.1: ESG migration matrix (Question 3.a)</caption><thead><tr><th>Rating</th><th>s<sub>k</sub></th><th>p<sub>k</sub></th><th colspan="4">Transition probability p<sub>j,k</sub></th><th>T(R<sub>k</sub>)</th></tr></thead><tbody><tr><td>D</td><td>-1.50</td><td>6.68%</td><td>92.96%</td><td>7.04%</td><td>0.00%</td><td>0.00%</td><td>7.04%</td></tr><tr><td>C</td><td>0.00</td><td>43.32%</td><td>1.31%</td><td>95.03%</td><td>3.66%</td><td>0.00%</td><td>4.97%</td></tr><tr><td>B</td><td>1.50</td><td>43.32%</td><td>0.00%</td><td>3.66%</td><td>95.03%</td><td>1.31%</td><td>4.97%</td></tr><tr><td>A</td><td></td><td>6.68%</td><td>0.00%</td><td>0.00%</td><td>7.04%</td><td>92.96%</td><td>7.04%</td></tr></tbody><tfoot><tr><td colspan="4">T(R<sub>1</sub>,...,R<sub>K</sub>)</td><td colspan="4">5.25%</td></tr></tfoot></table>

<table><caption>Table B.2: ESG migration matrix (Question 3.b)</caption><thead><tr><th>Rating</th><th>s<sub>k</sub></th><th>p<sub>k</sub></th><th colspan="4">Transition probability p<sub>j,k</sub></th><th>T(R<sub>k</sub>)</th></tr></thead><tbody><tr><td>D</td><td>-0.67</td><td>25.00%</td><td>95.15%</td><td>4.85%</td><td>0.00%</td><td>0.00%</td><td>4.85%</td></tr><tr><td>C</td><td>0.00</td><td>25.00%</td><td>5.27%</td><td>88.38%</td><td>6.35%</td><td>0.00%</td><td>11.62%</td></tr><tr><td>B</td><td>0.67</td><td>25.00%</td><td>0.00%</td><td>6.35%</td><td>88.38%</td><td>5.27%</td><td>11.62%</td></tr><tr><td>A</td><td></td><td>25.00%</td><td>0.00%</td><td>0.00%</td><td>4.85%</td><td>95.15%</td><td>4.85%</td></tr></tbody><tfoot><tr><td colspan="4">T(R<sub>1</sub>,...,R<sub>K</sub>)</td><td colspan="4">8.23%</td></tr></tfoot></table>

<table><caption>Table B.3: ESG migration matrix (Question 3.c)</caption><thead><tr><th>Rating</th><th>s<sub>k</sub></th><th>p<sub>k</sub></th><th colspan="4">Transition probability p<sub>j,k</sub></th><th>T(R<sub>k</sub>)</th></tr></thead><tbody><tr><td>D</td><td>-1.28</td><td>10.00%</td><td>93.54%</td><td>6.46%</td><td>0.00%</td><td>0.00%</td><td>6.46%</td></tr><tr><td>C</td><td>0.00</td><td>40.00%</td><td>1.89%</td><td>94.14%</td><td>3.97%</td><td>0.00%</td><td>5.86%</td></tr><tr><td>B</td><td>1.28</td><td>40.00%</td><td>0.00%</td><td>3.97%</td><td>94.14%</td><td>1.89%</td><td>5.86%</td></tr><tr><td>A</td><td></td><td>10.00%</td><td>0.00%</td><td>0.00%</td><td>6.46%</td><td>93.54%</td><td>6.46%</td></tr></tbody><tfoot><tr><td colspan="4">T(R<sub>1</sub>,...,R<sub>K</sub>)</td><td colspan="4">5.98%</td></tr></tfoot></table>

(h) The ESG migration matrices are given in Tables B.1, B.2 and B.3. We deduce that the ESG rating system defined in Question 3.a is the best rating system if we would like to reduce the monthly turnover of ESG ratings.

(i) The relationship between the parameter $\sigma$ and the turnover $T(R_1, \dots, R_K)$ is given in Figure B.6.

(j) The relationship between the number of notches $K$ and the turnover $T(R_1, \dots, R_K)$ is given in Figure B.7.

(k) An ESG rating system is mainly quantitative and highly depends on the mapping function. This is not the case of a credit rating system, which is mainly qualitative and discretionary. This explains that the turnover of an ESG rating system is higher than the turnover of a credit rating system. The stabilization of the ESG rating system implies to reduce the turnover $T(R_1, \dots, R_K)$, which depends on three factors: (1) the number of notches<fnref n="1" /> $K$; (2) the volatility $\sigma$ of score changes; (3) the design of the ESG rating system ($s_1, \dots, s_{K-1}$). The turnover $T(R_1, \dots, R_K)$ has a big impact on an ESG exclusion (or negative screening) policy, because it creates noisy short-term entry/exit positions that do not necessarily correspond to a decrease or increase of the long-term ESG risks.

<footnote n="1">¹This is why ESG rating systems have less notches than credit rating systems.</footnote>
<figure><img src="image_546.png" /><figcaption>Figure B.7: Relationship between $K$ and $T(\mathcal{R}_1, \dots, \mathcal{R}_K)$</figcaption></figure>

<figure><img src="image_547.png" /></figure>
B.1.3 Equity Portfolio optimization with ESG scores (Exercise 3.4.1)

1. We assume that the CAPM is valid.

(a) Using the CAPM, we have:

$$ \mu_i = r + \beta_i (\mu_m - r) $$

For instance, we have $\mu_1 = 1\% + 0.10 \times (6\% - 1\%) = 1.5\%$, $\mu_2 = 1\% + 0.30 \times 5\% = 2.5\%$. Finally, we obtain $\mu = (1.5\%, 2.5\%, 3.5\%, 5.5\%, 7.5\%, 11\%)$.

(b) We have:

$$ \Sigma = \sigma_m^2 \beta \beta^\top + D $$

where $D = \text{diag}(\tilde{\sigma}_1^2, \dots, \tilde{\sigma}_6^2)$. The numerical value of $\Sigma$ is:

$$ \Sigma = \begin{pmatrix} 293 & & & & & \\ 12 & 325 & & & & \\ 20 & 60 & 356 & & & \\ 36 & 108 & 180 & 424 & & \\ 52 & 156 & 260 & 468 & 797 & \\ 80 & 240 & 400 & 720 & 1040 & 1744 \end{pmatrix} \times 10^{-4} $$

The computation of the asset volatility is equal to $\sigma_i = \sqrt{\Sigma_{i,i}}$. We deduce that $\sigma = (17.12\%, 18.03\%, 18.87\%, 20.59\%, 28.23\%, 41.76\%)$. The formula of the correlation is $\rho_{i,j} = (\sigma_i \sigma_j)^{-1} \Sigma_{i,j}$. We obtain the following correlation matrix expressed in %:

$$ \mathbb{C} = \begin{pmatrix} 100.00 & & & & & \\ 3.89 & 100.00 & & & & \\ 6.19 & 17.64 & 100.00 & & & \\ 10.21 & 29.09 & 46.33 & 100.00 & & \\ 10.76 & 30.65 & 48.81 & 80.51 & 100.00 & \\ 11.19 & 31.88 & 50.76 & 83.73 & 88.21 & 100.00 \end{pmatrix} $$

(c) We have:

$$ w^* = \frac{\Sigma^{-1} (\mu - r \mathbf{1})}{\mathbf{1}^\top \Sigma^{-1} (\mu - r \mathbf{1})} = \begin{pmatrix} 0.94\% \\ 2.81\% \\ 5.28\% \\ 24.34\% \\ 29.06\% \\ 37.57\% \end{pmatrix} $$

We deduce that $\mu(w^*) = w^{*\top}\mu = 7.9201\%$ and $\sigma(w^*) = \sqrt{w^{*\top}\Sigma w^*} = 28.3487\%$. The Sharpe ratio is then equal to:

$$ \text{SR}(w^* | r) = \frac{7.9201\% - 1\%}{28.3487\%} = 0.2441 $$

Finally, the ESG score of the tangency portfolio is:

$$ \mathbf{S}(w^*) = \sum_{i=1}^{6} w_i^* \mathbf{S}_i = -2.0347 $$
(d) We have:

$$\beta_i(w^*) = \frac{\mathbf{e}_i^\top \Sigma w^*}{\sigma^2(w^*)}$$

We obtain:

$$\beta(w^*) = \begin{pmatrix} 0.0723 \\ 0.2168 \\ 0.3613 \\ 0.6503 \\ 0.9393 \\ 1.4451 \end{pmatrix}$$

The computation of the implied expected return $\tilde{\mu}_i = r + \beta_i(w^*) (\mu(w^*) - r)$ gives:

$$\tilde{\mu} = \begin{pmatrix} 1.50\% \\ 2.50\% \\ 3.50\% \\ 5.50\% \\ 7.50\% \\ 11.00\% \end{pmatrix}$$

(e) We notice that $\beta_i(w^*) \neq \beta_i(w_m)$ but the risk premia deduced from the tangency portfolio are the same as those computed with the market portfolio. In fact, the market portfolio cannot coincide with the tangency portfolio because of the cash component. Let us assume that the allocation of $w_m$ is equal to $\alpha$ of the tangency portfolio $w^*$ and $1 - \alpha$ of the risk-free asset. We deduce that:

$$\beta(w_m) = \frac{\Sigma w_m}{\sigma^2(w_m)} = \frac{\alpha \Sigma w^*}{\alpha^2 \sigma^2(w^*)} = \frac{1}{\alpha} \beta(w^*)$$

and:

$$\alpha = \frac{\beta_i(w^*)}{\beta_i(w_m)}$$

The computation gives $\alpha = 72.25\%$. The market portfolio $w_m$ is equal to 72.25% of the tangency portfolio $w^*$ and 27.75% of the risk-free asset. We have:

$$\begin{aligned} \mu(w_m) &= r + \alpha (\mu(w^*) - r) \\ &= 1\% + 72.25\% \times (7.9201\% - 1\%) \\ &= 6\% \end{aligned}$$

and:

$$\begin{aligned} \sigma(w_m) &= \alpha \sigma(w^*) \\ &= 72.25\% \times 28.3487\% \\ &= 20.48\% \end{aligned}$$

We deduce that:

$$\mathrm{SR}(w_m | r) = \frac{6\% - 1\%}{20.48\%} = 0.2441$$

Curiously, we do not obtain the true value of the Sharpe ratio:

$$\mathrm{SR}(w_m | r) = \frac{6\% - 1\%}{20\%} = 0.25$$
because the volatility $\sigma(w_m)$ is not equal to 20%. The reason is that the market portfolio computed with the tangency portfolio has an idiosyncratic risk. Indeed, we have:

$$ \sqrt{w_m^\top (\sigma_m^2 \beta \beta^\top) w^\top} = 20\% < \sigma(w_m) = 20.48\% $$

Therefore, $w_m$ is not fully diversified because the number of assets is small ($n = 6$).

2. We consider long-only portfolios and we also impose a minimum threshold $\mathbf{S}^*$ for the portfolio ESG score:

$$ \mathbf{S}(w) = \sum_{i=1}^{n} w_i \mathbf{S}_i \geq \mathbf{S}^* $$

(a) We have:

$$ \begin{array}{l} w^* = \arg \min \frac{1}{2} w^\top \Sigma w - \gamma w^\top \mu \\ \text{s.t.} \quad \begin{cases} \mathbf{1}_6^\top w = 1 \\ w^\top \mathbf{S} \geq \mathbf{S}^* \\ \mathbf{0}_6 \leq w \leq \mathbf{1}_6 \end{cases} \end{array} $$

<figure><img src="image_548.png" /><figcaption>Figure B.8: Impact of the minimum ESG score on the efficient frontier</figcaption></figure>

(b) We remind that the matrix form of the QP problem is:

$$ \begin{array}{l} w^* = \arg \min \frac{1}{2} w^\top Q w - w^\top R \\ \text{s.t.} \quad \begin{cases} Aw = B \\ Cw \leq D \\ w^- \leq w \leq w^+ \end{cases} \end{array} $$
We deduce that $Q = \Sigma$, $R = \gamma\mu$, $A = \mathbf{1}_6^\top$, $B = 1$, $C = -\mathbf{S}^\top$, $D = -\mathbf{S}^*$, $w^- = \mathbf{0}_6$ and $w^+ = \mathbf{1}_6$.

(c) To compute the efficient frontier, we consider several value of $\gamma \in [-1, 2]$. For each value of $\gamma$, we compute the optimal portfolio $w^*$ and deduce its expected return $\mu(w^*)$ and its volatility $\sigma(w^*)$. In Figure B.8, we compare the results for the three cases: $\mathbf{S}^* = -\infty$, $\mathbf{S}^* = 0$ and $\mathbf{S}^* = 0.5$. We notice that imposing a positive ESG score has little impact, especially for low-volatility portfolios when $\sigma(w^*) \le 16\%$. This is not the case when $\mathbf{S}^* = 0.5$.

(d) Let $w^*(\gamma)$ be the MVO portfolio when the risk tolerance is equal to $\gamma$. If we use a fine grid of $\gamma$ values, we can find the optimal value $\gamma^*$ by solving numerically the following optimization problem with the brute force algorithm:

$$ \gamma^* = \arg\max \frac{\mu(w^*(\gamma)) - r}{\sigma(w^*(\gamma))} \quad \text{for } \gamma \in [0, 2] $$

Then we deduce the tangency portfolio $w^* = w^*(\gamma^*)$. We also compute the following statistics: $\mu(w^*) = w^{*\top}\mu$, $\sigma(w^*) = \sqrt{w^{*\top}\Sigma w^*}$, $\text{SR}(w^*|r) = \sigma(w^*)^{-1}(\mu(w^*) - r)$ and $\mathbf{S}(w^*) = w^{*\top}\mathbf{S}$. Results are reported in Table B.4. In the case $\mathbf{S}^* = -\infty$, we retrieve the tangency portfolio, which has been found in Question 1.c. Moreover, we notice that the Sharpe ratio decreases when the minimum ESG score $\mathbf{S}^*$ increases.

<table><caption>Table B.4: Impact of the minimum ESG score on the efficient frontier</caption><thead><tr><th><b>S</b><sup>*</sup></th><th>-∞</th><th>0</th><th>0.5</th></tr></thead><tbody><tr><td>γ<sup>*</sup></td><td>1.1613</td><td>0.8500</td><td>0.8500</td></tr><tr><td rowspan="7">w<sup>*</sup> (in %)</td><td>0.9360</td><td>9.7432</td><td>9.1481</td></tr><tr><td>2.8079</td><td>16.3317</td><td>19.0206</td></tr><tr><td>5.2830</td><td>31.0176</td><td>40.3500</td></tr><tr><td>24.3441</td><td>5.1414</td><td>0.0000</td></tr><tr><td>29.0609</td><td>11.6028</td><td>3.8248</td></tr><tr><td>37.5681</td><td>26.1633</td><td>27.6565</td></tr><tr><td>μ(w<sup>*</sup>) (in %)</td><td>7.9201</td><td>5.6710</td><td>5.3541</td></tr><tr><td>σ(w<sup>*</sup>) (in %)</td><td>28.3487</td><td>19.8979</td><td>19.2112</td></tr><tr><td>SR(w<sup>*</sup>|r)</td><td>0.2441</td><td>0.2347</td><td>0.2266</td></tr><tr><td><b>S</b>(w<sup>*</sup>)</td><td>-2.0347</td><td>0.0000</td><td>0.5000</td></tr></tbody></table>

(e) We perform the same analysis as previously for several values $\mathbf{S}^* \in [-2.5, 2.5]$. Results are reported in Figure B.9. We verify that the Sharpe ratio is a decreasing function of $\mathbf{S}^*$.

(f) The market portfolio $w_m$ is then equal to:

$$ w_m = \begin{pmatrix} 9.15\% \\ 19.02\% \\ 40.35\% \\ 0.00\% \\ 3.82\% \\ 27.66\% \end{pmatrix} $$

We deduce that $\mu(w_m) = 5.3541\%$, $\sigma(w_m) = 19.2112\%$, $\text{SR}(w_m|r) = 0.2266$ and $\mathbf{S}(w_m) = 0.5$.
<figure><img src="image_549.png" /><figcaption>Figure B.9: Relationship between the minimum ESG score $\mathbf{S}^*$ and the Sharpe ratio SR ($w^* | r$) of the tangency portfolio</figcaption></figure>

(f).i We have:

$$ \beta_i(w_m) = \frac{\mathbf{e}_i^\top \Sigma w_m}{\sigma^2(w_m)} $$

and:

$$ \tilde{\mu}_i(w_m) = r + \beta_i(w_m) (\mu(w_m) - r) $$

We deduce that the alpha return is equal to:

$$ \begin{aligned} \alpha_i &= \mu_i - \tilde{\mu}_i(w_m) \\ &= (\mu_i - r) - \beta_i(w_m) (\mu(w_m) - r) \end{aligned} $$

Results are reported in Table B.5. We notice that $\alpha_i < 0$ for the first three assets and $\alpha_i > 0$ for the last three assets, implying that:

$$ \begin{cases} \mathbf{S}_i > 0 \Rightarrow \alpha_i < 0 \\ \mathbf{S}_i < 0 \Rightarrow \alpha_i > 0 \end{cases} $$

(f).ii We have:

$$ \beta(w_{\text{ew}} | w_m) = \frac{w_{\text{ew}}^{\top} \Sigma w_m}{\sigma^2(w_m)} = 0.9057 $$

and:

$$ \begin{aligned} \tilde{\mu}(w_{\text{ew}}) &= r + \beta(w_{\text{ew}} | w_m) (\mu(w_m) - r) \\ &= 1\% + 0.9057 \times (5.3541\% - 1\%) \\ &= 4.9435\% \end{aligned} $$