# -*- coding: utf-8 -*-
"""
This file estimates long run investment efficiencies related to increased certainty
or durability from the TPM proposal under a range of possible parameter values (i.e. monte-carlo)

Note that parameter distributions are calibrated so that the
mean values are approximately equal to the values outlined in the 
accompanying parameter assumption table in the document describing the methods 
used here. However the means produced by this file may not be exactly the same 
as in the parameter table.

Furthermore the quantity value listed in the method document 
(geometric mean MWh) is represented here by growth rates, 
for implementing the monte carlo. This requires a starting value 
which is the same starting values as used for the calculation of 
geometric mean quantity in the method document.

In the code below references to the parameters in the method document
are listed with e.g. 'Row 23', for row 23 in the table.

Calculations are also referenced to equations in the method document with e.g.
'Eq 17' for Equation 17 in the method document.

@author: theot
"""
#Libraries

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd
from scipy.stats import binom
from scipy.stats import beta
from random import random

### Efficient investment

#Unchanging parameters
T = 28 #Years 2022-2049
r = 0.06 #discount rate

P0 = 82 #Initial price, MBIE EDGS mixed renewables price index $/ MWh 2018
Q0 = 39572000 #Initial quantity  - MBIE, all consumption

seed_val = 999 #random seed value


# Parameters to vary
# Long run demand elasticity, Row 22
eta_alpha = 5 
eta_beta = 2
np.random.seed(seed=seed_val)
eta_rv = beta.rvs(eta_alpha, eta_beta, size=10000)
np.mean(eta_rv)
eta_p = eta_alpha, eta_beta #eta parameters
plt.figure()
plt.hist(eta_rv, 100, density=True, align='mid')
plt.xlabel('Absolute value of elasticity')
plt.ylabel('Frequency')
plt.title('Long run price elasticity of demand ($\\eta$), Beta distribution')
plt.text(0.3, 2, "$\\alpha$  = "+ str(eta_p[0]))
plt.text(0.3, 1.75, "$\\beta$ = " +str(eta_p[1]))



#Long run supply elasticity, Row 23
s_mu = 1 
s_sd = 0.25
np.random.seed(seed=seed_val)
s_rv = np.random.normal(s_mu, s_sd, 10000)
np.mean(s_rv)
s_p = s_mu, s_sd #supply elasticity parameters
plt.figure()
plt.hist(s_rv, 100, density=True, align='mid')
plt.xlabel('Supply elasticity')
plt.ylabel('Frequency')
plt.title('Long run supply elasticity, Normal distribution')
plt.text(1.5, 1.2, "$\\mu$  = "+ str(s_p[0]))
plt.text(1.5, 1.1, "$\\sigma$ = " +str(s_p[1]))

# Effect of investment on demand and supply, Row 28 
dQdI_mu = 0.05 
dQdI_sd = 0.01 
np.random.seed(seed=seed_val)
dQdI_rv = np.random.normal(dQdI_mu, dQdI_sd, 10000)
dQdI_p = dQdI_mu, dQdI_sd #dQdI parameters
plt.figure()
plt.hist(dQdI_rv, 100, density=True, align='mid')
plt.xlabel('Change in quantity given a change in investment')
plt.ylabel('Frequency')
plt.title('Effect of investment on demand and supply, Normal distribution')
plt.text(0.07, 30, "$\\mu$  = "+ str(dQdI_p[0]))
plt.text(0.07, 27.5, "$\\sigma$ = " +str(dQdI_p[1]))

#Effect of uncertainty on investment, Row 24
dIdU_mu = -0.087 
dIdU_sd = 0.01
np.random.seed(seed=seed_val) 
dIdU_rv = np.random.normal(dIdU_mu, dIdU_sd, 10000)
dIdU_p = dIdU_mu, dIdU_sd #dIdU parameters
plt.figure()
plt.hist(dIdU_rv, 100, density=True, align='mid')
plt.xlabel('Change in investment given a change in uncertainty')
plt.ylabel('Frequency')
plt.title('Effect of uncertainty on investment, Normal distribution')
plt.text(-0.065, 30, "$\\mu$  = "+ str(dIdU_p[0]))
plt.text(-0.065, 27.5, "$\\sigma$ = " +str(dIdU_p[1]))

#Baseline demand growth, Row 27
Qg_mu = 0.009 
Qg_sd = 0.005
np.random.seed(seed=seed_val) 
Qg_rv = np.random.normal(Qg_mu, Qg_sd, 10000)
Qg_p = Qg_mu, Qg_sd #Qg parameters
plt.figure()
plt.hist(Qg_rv, 100, density=True, align='mid')
plt.xlabel('Growth rate')
plt.ylabel('Frequency')
plt.title('Baseline demand growth, Normal distribution')
plt.text(0.02, 60, "$\\mu$  = "+ str(Qg_p[0]))
plt.text(0.02, 55, "$\\sigma$ = " +str(Qg_p[1]))

#Baseline price growth, Row 21 
Pg_mu = 0.00906 
Pg_sd = 0.0050 
np.random.seed(seed=seed_val)
Pg_rv = np.random.normal(Pg_mu, Pg_sd, 10000)
Pg_p = Pg_mu, Pg_sd #Pg parameters
plt.figure()
plt.hist(Pg_rv, 100, density=True, align='mid')
plt.xlabel('Growth rate')
plt.ylabel('Frequency')
plt.title('Baseline price growth, Normal distribution')
plt.text(0.02, 60, "$\\mu$  = "+ str(Pg_p[0]))
plt.text(0.02, 55, "$\\sigma$ = " +str(Pg_p[1]))

# Percentage change in frequency of events inducing uncertainty Row 26
dF_mu = -0.09 
dF_sd = 0.03 
np.random.seed(seed=seed_val)
dF_rv = np.random.normal(dF_mu, dF_sd, 10000)
dF_p = dF_mu, dF_sd #dF parameters
plt.figure()
plt.hist(dF_rv, 100, density=True, align='mid')
plt.xlabel('Percentage change in decimals')
plt.ylabel('Frequency')
plt.title('Change in frequency of events inducing uncertainty, Normal distribution')
plt.text(-0.04, 10, "$\\mu$  = "+ str(dF_p[0]))
plt.text(-0.04, 9, "$\\sigma$ = " +str(dF_p[1]))

 
# Cost benefit function
def f_benefit(eta_p = eta_p, s_p = s_p, dQdI_p = dQdI_p,
                       dIdU_p = dIdU_p, Qg_p = Qg_p, Pg_p = Pg_p,
                       dF_p = dF_p, Q0 = Q0, P0=P0):
    #Parameters
    eta = -1*beta.rvs(eta_p[0], eta_p[1])
    s =  np.random.normal(s_p[0], s_p[1])
    dQdI = np.random.normal(dQdI_p[0], dQdI_p[1]) # Eq 20
    dIdU = np.random.normal(dIdU_p[0], dIdU_p[1]) # Eq 20
    Qg =  np.random.normal(Qg_p[0], Qg_p[1])
    Pg =  np.random.normal(Pg_p[0], Pg_p[1])
    dF = np.random.normal(dF_p[0], dF_p[1])
    
    #Average demand and prices - baseline
    Q = (Q0*(Q0*(1+Qg)**T))**(0.5) #Row 21
    P = (P0*(P0*(1+Pg)**T))**(0.5) # Row 27

    #Supply and demand curve parameters
    Beta_s = s/(P/Q) # Eq 21, supply version
    Beta_d = eta/(P/Q) # Eq 21
    Q_P0 = Q+(-1*Beta_d*P) # Eq 19
    P_Q0 = Q_P0/(-1*Beta_d) # Eq 18
    dU = (dIdU*dQdI)/(P/Q) # Eq 20
    dPdU = (dU*2)/(Beta_s-Beta_d) # Eq 15
    dP = dF*dPdU
    dQ = dF*100*dU
    dTS_pa = 0.5*(((P*(1+dP)*(Q+dQ))-P*Q)+(((P_Q0-P*(1+dP))*(Q+dQ))-(P_Q0-P)*Q)) #Eq 17
    dTS = ((dTS_pa/r)*(1-(1/(1+r)**T)))/(1+r)**2 # Note additional deflator term accounts for time between 2019 and 2022
    
    return dTS

# Simulation - random
sims = 50000
# Initialise
dTS = np.zeros(sims)
#Evaluate
for i in range(sims):
    dTS[i] = f_benefit()
    
matplotlib.pyplot.close("all")
plt.figure()
plt.hist(dTS/1000000,100, density=True, align='mid') 
plt.xlabel('Change in total surplus ($ millions present value)')
plt.ylabel('Frequency')
plt.title('Durability improvement, monte carlo results')
plt.text(80, 0.02, 'Mean = %s'%(round(np.mean(dTS)/1000000,2)))

np.percentile(dTS,5)
np.percentile(dTS,50)
np.percentile(dTS,95)
np.mean(dTS)

ptiles = [[np.percentile(dTS,5),np.percentile(dTS,50),np.mean(dTS),np.percentile(dTS,95)]]
ptiles = pd.DataFrame(ptiles)
ptiles.columns = ["5th percentile","50th percentile","Mean","95th percentile"]
ptiles.to_csv('Durability_results_final.csv')