.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_data_analysis/distribution_fitting/plot_estimate_gpd_dowjones.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_data_analysis_distribution_fitting_plot_estimate_gpd_dowjones.py: Estimate a GPD on the Dow Jones Index data ========================================== .. GENERATED FROM PYTHON SOURCE LINES 6-9 In this example, we illustrate various techniques of extreme value modeling applied to the 5-year series of daily Dow Jones Index closing prices. Readers should refer to [coles2001]_ example 1.8 to get more details. .. GENERATED FROM PYTHON SOURCE LINES 9-16 .. code-block:: Python import math as m import openturns as ot import openturns.experimental as otexp import openturns.viewer as otv from openturns.usecases import coles import pandas as pd .. GENERATED FROM PYTHON SOURCE LINES 17-19 First, we load the Dow Jones dataset and plot it through time. We can see that the process is non-stationary. .. GENERATED FROM PYTHON SOURCE LINES 19-34 .. code-block:: Python full = pd.read_csv(coles.Coles().dowjones, index_col=0, parse_dates=True) print(full[:10]) graph = ot.Graph( "Daily closing prices of the Dow Jones Index", "Day index", "Index", True, "" ) # Care: to get the real period range, multiply by a factor to account for non-working days missing values! size = len(full) days = ot.Sample([[i] for i in range(size)]) dataDowJones = ot.Sample.BuildFromDataFrame(full) curve = ot.Curve(days, dataDowJones) curve.setColor("red") graph.add(curve) graph.setIntegerXTick(True) view = otv.View(graph) .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_dowjones_001.png :alt: Daily closing prices of the Dow Jones Index :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_dowjones_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Index Date 1995-09-11 02:00:00 4704.94 1995-09-12 02:00:00 4747.21 1995-09-13 02:00:00 4765.52 1995-09-14 02:00:00 4801.80 1995-09-15 02:00:00 4797.57 1995-09-18 02:00:00 4780.41 1995-09-19 02:00:00 4767.04 1995-09-20 02:00:00 4792.69 1995-09-21 02:00:00 4767.40 1995-09-22 02:00:00 4764.15 .. GENERATED FROM PYTHON SOURCE LINES 35-45 In that example, the time dependence can not be explained by trends or seasonal cycles. Many empirical studies have advised to consider the logarithms of ratios of successive observations to get an approximation to stationarity. We apply that transformation: .. math:: \tilde{X}_i = \log X_i - \log X_{i-1}. The resulting time series appears to be reasonably close to stationarity. .. GENERATED FROM PYTHON SOURCE LINES 45-59 .. code-block:: Python transfDataDJ = ot.Sample( [ [m.log(dataDowJones[i, 0]) - m.log(dataDowJones[i - 1, 0])] for i in range(1, size) ] ) curve = ot.Curve(days[:-1], transfDataDJ) graph = ot.Graph( "Log-daily returns of the Dow Jones Index", "Day index", "Index", True, "" ) graph.add(curve) view = otv.View(graph) .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_dowjones_002.png :alt: Log-daily returns of the Dow Jones Index :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_dowjones_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 60-65 For convenience of presentation, we rescale the data: .. math:: \tilde{X}_i \rightarrow 100\tilde{X}_i. .. GENERATED FROM PYTHON SOURCE LINES 65-68 .. code-block:: Python scalTransfDataDJ = transfDataDJ * 100.0 size = len(scalTransfDataDJ) .. GENERATED FROM PYTHON SOURCE LINES 69-75 In order to select a threshold upon which the GPD model will be fitted, we draw the mean residual life plot with approximate :math:`95\%` confidence interval. The curve is initially linear and shows significant curvature for :math:`u \in [-1, 2]`. Then for :math:`u \geq 2`, the curve is considered as reasonably linear when judged to confidence intervals. Hence, we choose the threshold :math:`u=2`. There are 37 exceedances of :math:`u`. .. GENERATED FROM PYTHON SOURCE LINES 75-80 .. code-block:: Python factory = ot.GeneralizedParetoFactory() graph = factory.drawMeanResidualLife(scalTransfDataDJ) view = otv.View(graph) u = 2.0 .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_dowjones_003.png :alt: Mean residual life plot :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_dowjones_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 81-88 **Stationary GPD modeling assuming independence in data** We first assume that the dependence between the transformed data is negligible, so we first consider the data as independent observations over the observation period. We estimate the parameters of the GPD distribution modeling the excesses above :math:`u=2` by maximizing the log-likelihood of the excesses. .. GENERATED FROM PYTHON SOURCE LINES 88-90 .. code-block:: Python result_LL = factory.buildMethodOfLikelihoodMaximizationEstimator(scalTransfDataDJ, u) .. GENERATED FROM PYTHON SOURCE LINES 91-92 We get the fitted GPD and the estimated parameters :math:`(\hat{\sigma}, \hat{\xi})`. .. GENERATED FROM PYTHON SOURCE LINES 92-98 .. code-block:: Python fitted_GPD = result_LL.getDistribution() desc = fitted_GPD.getParameterDescription() param = fitted_GPD.getParameter() print(", ".join([f"{p}: {value:.3f}" for p, value in zip(desc, param)])) print("log-likelihood = ", result_LL.getLogLikelihood()) .. rst-class:: sphx-glr-script-out .. code-block:: none sigma: 0.626, xi: 0.066, u: 2.000 log-likelihood = -22.103513585430218 .. GENERATED FROM PYTHON SOURCE LINES 99-106 We get the asymptotic distribution of the estimator :math:`(\hat{\sigma}, \hat{\xi})`. The threshold :math:`u` has not been estimated to ensure the regularity of the model and then the asymptotic properties of the maximum likelihood estimators. This is the reason why it appears as a Dirac distribution centered on the chosen threshold. In that case, the asymptotic distribution of :math:`(\hat{\sigma}, \hat{\xi})` is normal. .. GENERATED FROM PYTHON SOURCE LINES 106-110 .. code-block:: Python parameterEstimate = result_LL.getParameterDistribution() print("Asymptotic distribution of the estimator : ") print(parameterEstimate) .. rst-class:: sphx-glr-script-out .. code-block:: none Asymptotic distribution of the estimator : BlockIndependentDistribution(Normal(mu = [0.625758,0.0661835], sigma = [0.1838,0.193486], R = [[ 1 -0.795905 ] [ -0.795905 1 ]]), Dirac(point = [2])) .. GENERATED FROM PYTHON SOURCE LINES 111-113 We get the covariance matrix and the standard deviation of :math:`(\hat{\sigma}, \hat{\xi}, u)` where :math:`u` is deterministic. .. GENERATED FROM PYTHON SOURCE LINES 113-116 .. code-block:: Python print("Cov matrix = \n", parameterEstimate.getCovariance()) print("Standard dev = ", parameterEstimate.getStandardDeviation()) .. rst-class:: sphx-glr-script-out .. code-block:: none Cov matrix = [[ 0.0337823 -0.0283045 0 ] [ -0.0283045 0.0374369 0 ] [ 0 0 0 ]] Standard dev = [0.1838,0.193486,0] .. GENERATED FROM PYTHON SOURCE LINES 117-119 We get the marginal confidence intervals of order 0.95 of :math:`(\hat{\sigma}, \hat{\xi})`. .. GENERATED FROM PYTHON SOURCE LINES 119-124 .. code-block:: Python order = 0.95 for i in range(2): # exclude u parameter (fixed) ci = parameterEstimate.getMarginal(i).computeBilateralConfidenceInterval(order) print(desc[i] + ":", ci) .. rst-class:: sphx-glr-script-out .. code-block:: none sigma: [0.265518, 0.985999] xi: [-0.313043, 0.44541] .. GENERATED FROM PYTHON SOURCE LINES 125-137 At last, we can check the quality of the inference thanks to the 4 usual diagnostic plots: - the probability-probability pot, - the quantile-quantile pot, - the return level plot, - the data histogram and the density of the fitted model. We conclude that the goodness-of-fit in the quantile plots seems unconvincing, even if the other plots appear to be reasonable. This is due to the fact that the excesses can not be considered as independent: the transformed series :math:`\tilde{X}_i` has a rich structure of temporal dependence. .. GENERATED FROM PYTHON SOURCE LINES 137-141 .. code-block:: Python validation = otexp.GeneralizedParetoValidation(result_LL, scalTransfDataDJ) graph = validation.drawDiagnosticPlot() view = otv.View(graph) .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_dowjones_004.png :alt: , Sample versus model PP-plot, Sample versus model QQ-plot, Return level plot, Density :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_dowjones_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 142-160 **Stationary GPD modeling taking into account the dependence in data** We illustrate the fact that the excesses of the transformed series happen in groups. Hence we use the declustering method which filters the dependent observations exceeding a given threshold to obtain a set of threshold excesses that can be assumed as independent. Consecutive exceedances of :math:`u` belong to the same cluster. Two distinct clusters are separated by :math:`r` consecutive observations under the threshold. Within each cluster, we select the maximum value that will be used to infer the GPD distribution. The cluster maxima are assumed to be independent. On the graph, we show the clusters over the threshold :math:`u=2` and all the maxima selected within each cluster. It is possible to extract the data belonging to the same cluster and the cluster maximum series. We denote by :math:`n_c` the number of clusters and by :math:`n_u` the number of exceedances above :math:`u`. .. GENERATED FROM PYTHON SOURCE LINES 160-172 .. code-block:: Python part = otexp.SamplePartition(scalTransfDataDJ) r = 3 peaks, clusters = part.getPeakOverThreshold(u, r) nc = len(peaks) nu = sum([1 if scalTransfDataDJ[i, 0] > u else 0 for i in range(size)]) print(f"nc={nc} nu={u} theta={nc / nu:.3f}") graph = clusters.draw(u) graph.setTitle( "Threshold exceedances and clusters by transformed Dow Jones Index series" ) view = otv.View(graph) .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_dowjones_005.png :alt: Threshold exceedances and clusters by transformed Dow Jones Index series :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_dowjones_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none nc=32 nu=2.0 theta=0.865 .. GENERATED FROM PYTHON SOURCE LINES 173-175 We estimate a stationary GPD on the clusters maxima which are independent with the :math:`95\%` confidence interval of each parameter. .. GENERATED FROM PYTHON SOURCE LINES 175-183 .. code-block:: Python result_LL = factory.buildMethodOfLikelihoodMaximizationEstimator(peaks, u) sigma, xi, _ = result_LL.getParameterDistribution().getMean() sigma_stddev, xi_stddev, _ = result_LL.getParameterDistribution().getStandardDeviation() print( f"u={u} r={r} nc={nc} sigma={sigma:.2f} ({sigma_stddev:.2f}) xi={xi:.2f} ({xi_stddev:.2f})", end=" ", ) .. rst-class:: sphx-glr-script-out .. code-block:: none u=2.0 r=3 nc=32 sigma=0.54 (0.18) xi=0.27 (0.28) .. GENERATED FROM PYTHON SOURCE LINES 184-191 We evaluate the :math:`T=100`-year return level which corresponds to the :math:`m`-observation return level, where :math:`m = T*n_y` with :math:`n_y` the number of observations per year. Here, we have daily observations, hence :math:`n_y = 365`. To calculate it, we evaluate the extremal index :math:`\theta` which is the inverse of the mean length of the clusters, estimated by the ratio between the number of clusters and the number of exceedances of :math:`u`. .. GENERATED FROM PYTHON SOURCE LINES 191-197 .. code-block:: Python theta = nc / nu ny = 365 T = 100 xm_100 = factory.buildReturnLevelEstimator(result_LL, scalTransfDataDJ, T * ny, theta) print(f"x100={xm_100.getMean()} ({xm_100.getStandardDeviation()}) theta={theta:.3f}") .. rst-class:: sphx-glr-script-out .. code-block:: none x100=[12.5097] ([10.6691]) theta=0.865 .. GENERATED FROM PYTHON SOURCE LINES 198-203 We plot the return level for the new fitted model that takes into account dependence between excesses. We can see the fitted model works well. However, the large return level confidence intervals obtained for extreme return levels makes it difficult to make reliable predictions with any degree of certainty. .. GENERATED FROM PYTHON SOURCE LINES 203-209 .. code-block:: Python validation = otexp.GeneralizedParetoValidation(result_LL, peaks) grid = validation.drawDiagnosticPlot() rlPlot = grid.getGraph(1, 0) rlPlot.setTitle(rlPlot.getTitle() + f" (u={u} r={r})") view = otv.View(rlPlot) .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_dowjones_006.png :alt: Return level plot (u=2.0 r=3) :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_dowjones_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 210-211 We plot the whole series of validation graphs of the new fitted model. .. GENERATED FROM PYTHON SOURCE LINES 211-213 .. code-block:: Python view = otv.View(grid) .. image-sg:: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_dowjones_007.png :alt: , Sample versus model PP-plot, Sample versus model QQ-plot, Return level plot, Density :srcset: /auto_data_analysis/distribution_fitting/images/sphx_glr_plot_estimate_gpd_dowjones_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 214-215 .. code-block:: Python otv.View.ShowAll() .. _sphx_glr_download_auto_data_analysis_distribution_fitting_plot_estimate_gpd_dowjones.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_estimate_gpd_dowjones.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_estimate_gpd_dowjones.py `