code: #!/usr/bin/env python3
  # -*- coding: utf-8 -*-
  # ---
  # jupyter:
  #   jupytext:
  #     text_representation:
  #       extension: .py
  #       format_name: light
  #       format_version: '1.4'
  #       jupytext_version: 1.2.0
  #   kernelspec:
  #     display_name: Python 3
  #     language: python
  #     name: python3
  # ---
  
  # # s_checklist_scenariobased_step07 [<img src="https://www.arpm.co/lab/icons/icon_permalink.png" width=30 height=30 style="display: inline;">](https://www.arpm.co/lab/redirect.php?code=s_checklist_scenariobased_step07&codeLang=Python)
  # For details, see [here](https://www.arpm.co/lab/redirect.php?permalink=ex-vue-7).
  
  # +
  import numpy as np
  import pandas as pd
  
  from arpym.portfolio import spectral_index
  from arpym.statistics import meancov_sp, quantile_sp
  
  # -
  
  # ## [Input parameters](https://www.arpm.co/lab/redirect.php?permalink=s_checklist_scenariobased_step07-parameters)
  
  # +
  # indicates which projection to continue from
  # True: use copula-marginal projections
  # False: use historical projections
  copula_marginal = True
  
  lam = 3e-7  # parameter of exponential utility function
  c_quantile = 0.95  # confidence level for the quantile satisfaction measure
  c_es = 0.95  # confidence level for the negative expected shortfall
  # -
  
  # ## [Step 0](https://www.arpm.co/lab/redirect.php?permalink=s_checklist_scenariobased_step07-implementation-step00): Load data
  
  # +
  path = '../../../databases/temporary-databases/'
  
  if copula_marginal:
      # Projection
      db_projection_tools = pd.read_csv(path + 'db_projection_tools.csv')
      j_ = int(db_projection_tools['j_'][0])
  
      db_scenprob = pd.read_csv(path + 'db_scenario_probs.csv')
      p = db_scenprob['p'].values
  
      # Pricing
      db_pricing = pd.read_csv(path + 'db_pricing.csv')
      pi_tnow_thor = db_pricing.values
  
      # Aggregation
      db_exante_perf = pd.read_csv(path + 'db_exante_perf.csv')
      y_h = db_exante_perf.values.squeeze()
  
  else:
      # Projection
      db_projection_tools = pd.read_csv(path + 'db_projection_bootstrap_tools.csv')
      j_ = int(db_projection_tools['j_'][0])
  
      db_scenprob = pd.read_csv(path + 'db_scenario_probs_bootstrap.csv')
      p = db_scenprob['p'].values
  
      # Pricing
      db_pricing = pd.read_csv(path + 'db_pricing_historical.csv')
      pi_tnow_thor = db_pricing.values
  
      # Aggregation
      db_exante_perf = pd.read_csv(path + 'db_exante_perf_historical.csv')
      y_h = db_exante_perf.values.squeeze()
  
  db_holdings = pd.read_csv(path + 'db_holdings.csv')
  h = np.squeeze(db_holdings.values)
  # -
  
  # ## [Step 1](https://www.arpm.co/lab/redirect.php?permalink=s_checklist_scenariobased_step07-implementation-step01): Calculate certainty equivalent satisfaction measure
  
  # +
  # expected utility
  expected_utility = p @ (-np.exp(-lam * y_h))  # expected utility computation
  
  # certainty equivalent satisfaction measure
  cert_eq_yh = -(1 / lam) * np.log(-expected_utility)
  # -
  
  # ## [Step 2](https://www.arpm.co/lab/redirect.php?permalink=s_checklist_scenariobased_step07-implementation-step02): Quantile satisfaction measure
  
  # quantile
  q_yh = quantile_sp(1 - c_quantile, y_h, p, method='kernel_smoothing')
  
  
  # ## [Step 3](https://www.arpm.co/lab/redirect.php?permalink=s_checklist_scenariobased_step07-implementation-step03): Expected shortfall satisfaction measure
  
  # +
  # indicator function
  def indicator(x):
      return (0 <= x and x <= 1 - c_es)
  
  
  # spectrum function
  def spectr_es(x):
      return (1 / (1 - c_es)) * indicator(x)
  
  
  # negative expected shortfall
  es_yh, _ = spectral_index(spectr_es, pi_tnow_thor,
                            p, h)
  # -
  
  # ## [Step 4](https://www.arpm.co/lab/redirect.php?permalink=s_checklist_scenariobased_step07-implementation-step04): Expectation and variance satisfaction measures
  
  # expectation satisfaction measure
  mean_yh, var_yh = meancov_sp(y_h, p)
  # opposite of variance is satisfaction measure
  neg_var_yh = -var_yh
  
  # ## [Step 5](https://www.arpm.co/lab/redirect.php?permalink=s_checklist_scenariobased_step07-implementation-step05): Save database
  
  # +
  out = pd.DataFrame({'cert_eq_yh': pd.Series(cert_eq_yh),
                      'q_yh': pd.Series(q_yh),
                      'es_yh': pd.Series(es_yh),
                      'mean_yh': pd.Series(mean_yh),
                      'neg_var_yh': pd.Series(neg_var_yh),
                      'c_es': pd.Series(c_es),
                      'c_quantile': pd.Series(c_quantile)})
  if copula_marginal:
      out.to_csv(path + 'db_quantile_and_satis.csv',
                 index=False)
  else:
      out.to_csv(path + 'db_quantile_and_satis_historical.csv',
                 index=False)
  
  del out
   import matplotlib.pyplot as plt
  from numpy import array, ones, zeros, cov
  
  plt.style.use('seaborn')
  from fglasso import glasso
  
  
  def GraphicalLasso(pop, lam, initStruct=None, approximate=0, warmInit=0, verbose=0, penalDiag=1, tolThreshold=1e-4, maxIter=1e4, w=None, theta=None):
      # [w, theta, iter, avgTol, hasError] = GraphicalLasso(pop, lam,
      # initStruct, approximate, warmInit, verbose, penalDiag, tolThreshold,
      # maxIter, w, theta)
      #
      # Computes a regularized estimate of covariance matrix and its inverse
      # Inputs:
      #  - pop: the set of samples to be used for covariance estimation in an NxP
      #  matrix where N is the number of samples and P is the number of variables
      #  - lam: the regularization penalty. Can be a single number or a matrix
      #  of penalization values for each entry in the covariance matrix
      #  - initStruct(o@): a matrix of size PxP, where zero entries will force
      #  the corresponding entries in the inverse covariance matrix to be zero.
      #  - approximate([o]): a flag indicating whether to use approximate estimation
      #  (Meinhausen-Buhlmann approximation)
      #  - warmInit[o]: a flag indicating whether the estimation will start from
      #  given initial values of w and theta
      #  - verbose([o]): a flag indicating whether to output algorithm process
      #  - penalDiag[o]: a flag indicating whether to penalize diagonal elements of
      #  the covariance matrix
      #  - tolThreshold[o]: the amount of tolerance acceptable for covariance matrix
      #  elements before terminating the algorithm
      #  - maxIter([o]): maximum number of iteration to perform in the algorithm
      #  - w[o]: the initial value of covariance matrix used for warm initialization
      #  - theta([o]): the initial value of inverse covariance matrix used for warm
      #  initialization
      # @: o indicates optional arguments
      # Outputs:
      #  - w: the estimated covariance matrix
      #  - theta: the estimated inverse covariance matrix
      #  - iter: actual number of iterations performed in the algorithm
      #  - avgTol: average tolerance of covariance matrix entries before
      #  terminating the algorithm
      #  - hasError: a flag indicating whether the algorithm terminated
      #  erroneously or not
      #
      # Code by: Hossein Karshenas (hkarshenas@fi.upm.es)
      # Date: 10 Feb 2011
  
      numVars = pop.shape[1]
      if isinstance(lam,float):
          lam = array([[lam]])
      m, n = lam.shape
      if m != n:
          raise ValueError('Regularization coefficients matrix should be symmetric matrix')
      elif m > 1 and m != numVars:
          raise ValueError('Regularization coefficients matrix should have a size equal to the number of variables')
      if m == 1:
          lam = lam * ones((numVars,numVars))
  
      if initStruct is not None:
          initStruct = 1 - initStruct
          initStruct = 10e9*initStruct
          lam = lam + initStruct
  
      if w is None:
          if warmInit is False:
              raise ValueError('In warm initialization mode starting values for the covariance and precision matrices should be determined')
          else:
              w = zeros((numVars,numVars),order='F')
      if theta is None:
          if warmInit is False:
              raise ValueError('In warm initialization mode starting values for the precision matrix should be determined')
          else:
              theta = zeros((numVars,numVars),order='F')
  
      niter = 0
      jerr = 0.0
      # glasso(cov(pop.T), lam, 1, 1, 1, penalDiag, tolThreshold, maxIter, w, theta, niter, jerr, 1)
      glasso(cov(pop.T), lam, approximate, 0, verbose, penalDiag, tolThreshold, maxIter, w, theta, niter, jerr, numVars)
      # w, theta, iter, avgTol, hasError = glasso(cov(pop.T), lam, approximate, 0, verbose, penalDiag, tolThreshold, maxIter, w, theta, 0, 0, numVars)
  
      if False:
          raise Warning('The execution of the algorithm caused errors')
      return w, theta, niter, jerr
  
   #!/usr/bin/env python3
  # -*- coding: utf-8 -*-
  # ---
  # jupyter:
  #   jupytext:
  #     text_representation:
  #       extension: .py
  #       format_name: light
  #       format_version: '1.4'
  #       jupytext_version: 1.1.4
  #   kernelspec:
  #     display_name: Python 3
  #     language: python
  #     name: python3
  # ---
  
  # # S_PricingEquityTaylor [<img src="https://www.arpm.co/lab/icons/icon_permalink.png" width=30 height=30 style="display: inline;">](https://www.arpm.co/lab/redirect.php?code=S_PricingEquityTaylor&codeLang=Python)
  # For details, see [here](https://www.arpm.co/lab/redirect.php?permalink=eb-taylor-equity-pl).
  
  # ## Prepare the environment
  
  # +
  import os
  import os.path as path
  import sys
  
  sys.path.append(path.abspath('../../functions-legacy'))
  from collections import namedtuple
  
  from numpy import arange, array, ones, zeros, sort, where, diff, linspace, round, log, exp, sqrt
  from numpy import sum as npsum, min as npmin, max as npmax
  
  from scipy.stats import norm
  from scipy.io import loadmat
  
  import matplotlib.pyplot as plt
  from matplotlib.pyplot import figure, plot, legend, xlim, ylabel, \
      xlabel, title, xticks, yticks
  
  plt.style.use('seaborn')
  
  from CONFIG import GLOBAL_DB, TEMPORARY_DB
  from ARPM_utils import save_plot, struct_to_dict
  from FPmeancov import FPmeancov
  from HistogramFP import HistogramFP
  from SaddlePointQuadN import SaddlePointQuadN
  from SimulateBrownMot import SimulateBrownMot
  # -
  
  # ## run S_PricingEquityProfitAndLoss
  
  from S_PricingEquityProfitAndLoss import *
  
  # ## Compute the first order approximation of the equity P&L, which has normal distribution,
  # ## and the second order approximation of the equity P&L, that has generalized non-central gamma distribution,
  # ## at the selected horizon (120 days). Use function SaddlePointQuadN to compute the cumulative distribution
  # ## function of the generalized non-central gamma distribution.
  
  # +
  n_ = 500
  hor_sel = 120  # selected horizon for the plot (120 days)
  i = where(horiz_u == hor_sel)[0][-1]
  x_hor = zeros((n_, i+1))
  Taylor_first = zeros((n_, i+1))
  cdf_QuadN = zeros((n_, i+1))
  Taylor_second = zeros((n_, i+1))
  
  x_hor[:,i] = linspace(Mu_PL[0, i] - 10*Sigma_PL[0, i], Mu_PL[0, i] + 10 *Sigma_PL[0, i], n_)
  # first order approximation (normal pdf)
  Taylor_first[:,i] = norm.pdf(x_hor[:,i], exp(x[0,-1])*mu*horiz_u[i + 1],exp(x[0,-1])*sig*sqrt(horiz_u[i + 1]))
  # second order approximation (QuadN pdf)
  b, c, mu2, sigma2 = array([[exp(x[0,-1])]]), array([[exp(x[0,-1])*0.5]]), mu*horiz_u[i + 1], sig**2*horiz_u[i + 1]
  _, Taylor_second[:,i] = SaddlePointQuadN(x_hor[:,[i]].T, 0, b, c, mu2, sigma2)  # QuadN cumulative density function
  # Taylor_second(:,i) = diff(cdf_QuadN(:,i))/diff(x_hor((:,i)))
  # -
  
  # ## Plot a few (say 15) simulated paths of the equity P&L up to the selected horizon (120 days),
  # ## along with the first order approximation, the second order approximation and the analytical
  # ## distribution of the equity P&L. Furthermore, show the mean and the standard deviation of
  # ## the analytical distribution.
  
  # +
  lgrey = [0.8, 0.8, 0.8]  # light grey
  dgrey = [0.4, 0.4, 0.4]  # dark grey
  lblue = [0.27, 0.4, 0.9]  # light blu
  orange = [0.94, 0.35, 0]  # orange
  j_sel = 15  # selected MC simulations
  
  figure()
  # simulated path, mean and standard deviation
  plot(horiz_u[:i], PL[:j_sel,:i].T, color=lgrey)
  plt.xticks(arange(0,t_end+20,20))
  xlim([npmin(horiz_u), npmax(horiz_u)+1])
  l1 = plot(horiz_u[:i], Mu_PL[0,:i], color='g', label='mean')
  l2 = plot(horiz_u[:i], Mu_PL[0,:i] + Sigma_PL[0,:i], color='r', label=' + / - st.deviation')
  plot(horiz_u[:i], Mu_PL[0,:i] - Sigma_PL[0,:i], color='r')
  
  # analytical pdf
  flex_probs_scenarios = ones((1, j_)) / j_
  option = namedtuple('option','n_bins')
  option = namedtuple('option', 'n_bins')
  option.n_bins = round(10 * log(j_))
  y_hist, x_hist = HistogramFP(PL[:,[i]].T, flex_probs_scenarios, option)
  scale = 0.15 * Sigma_PL[0, i] / npmax(y_hist)
  y_hist = y_hist * scale
  shift_y_hist = horiz_u[i] + y_hist
  emp_pdf = plt.barh(x_hist[:-1], shift_y_hist[0]-npmin(shift_y_hist[0]), height=x_hist[1]-x_hist[0],
           left=npmin(shift_y_hist[0]), facecolor=lgrey, edgecolor=lgrey,label='horizon pdf')  # empirical pdf
  plot(shift_y_hist[0], x_hist[:-1], color=dgrey)  # border
  
  # first order approximation
  Taylor_first[:,i] = Taylor_first[:,i]*scale
  shift_T_first = horiz_u[i] + Taylor_first[:,i]
  l3 = plot(shift_T_first, x_hor[:,i], color=orange, label='first order approx')
  
  # second order approximation
  Taylor_second[:,i] = Taylor_second[:,i]*scale
  shift_T_second = horiz_u[i] + Taylor_second[:,i]
  l4 = plot(shift_T_second, x_hor[:,i], color=lblue, label='second order approx')
  
  legend()
  xlabel('time (days)')
  ylabel('P&L')
  title('P&L equity Taylor approximation');
  # save_plot(ax=plt.gca(), extension='png', scriptname=os.path.basename('.')[:-3], count=plt.get_fignums()[-1])
  
apis: pandas.Series numpy.squeeze pandas.read_csv numpy.exp numpy.log
apis: numpy.ones matplotlib.pyplot.style.use numpy.zeros numpy.array numpy.cov
apis: numpy.ones matplotlib.pyplot.style.use matplotlib.pyplot.legend numpy.zeros matplotlib.pyplot.figure numpy.exp matplotlib.pyplot.title numpy.arange numpy.max matplotlib.pyplot.ylabel numpy.min numpy.log numpy.sqrt matplotlib.pyplot.plot numpy.where numpy.linspace matplotlib.pyplot.xlabel