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

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
from scipy.stats import binom
from scipy.stats import beta
from random import random
import pandas as pd
### Efficient investment

#Unchanging parameters
T = 29 #Years
t0 = 2019 #Base year for discounting purpose
d = 1/1.06 #Discount rate
d_all = [1 * d ** (t - 1) for t in range(1, T + 3)] # Discount rates 2019-2049
dt=d_all[3:] #discount rates 2022-2049

I_0 = [27685076,80734170,43719050,67809249,55754905,15367492,15369216,
              17166939,45312958,164259283,81091654,171015884,140691756,
              89203879,137003955] # Forecast investment, E&D, 2021-2035
g_0 = 0.00884 #Forecast demand growth
Q_0=7376 # Initial year post diversity national peak demand
Q = np.full(T,Q_0) # Forecast post-diversity national peak demand 
for i in range(1,T):
    Q[i]=Q[i-1]*(1+g_0)
dQt = np.diff(Q) #Change in peak demand
c = 280617 # average investment cost per MW 2021-2035
I=np.concatenate((I_0,(np.multiply(dQt[len(I_0):,],c)))) # investment
ct = np.divide(I,dQt) #Unit investment costs over time
Qt = Q[1:]
seed_val = 999 #random seed value

# Parameters to vary
#Typical MW in area of benefit, share of peak, Row 13
QA_alpha = 2 
QA_beta = 4 
np.random.seed(seed=seed_val)
QA_rv = beta.rvs(QA_alpha, QA_beta, size=10000)
np.mean(QA_rv)
QA_p = QA_alpha, QA_beta # QA parameters

plt.figure()
plt.hist(QA_rv, 100, density=True, align='mid')
plt.xlabel('Share of peak MW')
plt.ylabel('Frequency')
plt.title('Typical MW in area of benefit ($Q_A$), Beta distribution')
plt.text(0.7, 1.5, "$\\alpha$  = "+ str(QA_p[0]))
plt.text(0.7, 1.25, "$\\beta$ = " +str(QA_p[1]))

# Transmission costs, share of prices
sT_mu = 0.0084 
sT_sd = 0.0005
np.random.seed(seed=seed_val) 
sT_rv = np.random.normal(sT_mu, sT_sd, 10000)
plt.figure()
plt.hist(sT_rv, 100, density=True, align='mid')
sT_p = sT_mu, sT_sd #sT parameters
plt.figure()
plt.hist(sT_rv, 100, density=True, align='mid')
plt.xlabel('Share of prices')
plt.ylabel('Frequency')
plt.title('Incremental transmission cost share of prices ($s_T$), Normal distribution')
plt.text(0.009, 600, "$\\mu$  = "+ str(sT_p[0]))
plt.text(0.009, 500, "$\\sigma$ = " +str(sT_p[1]))


# Share of investment for reliability, Row 2
sr_alpha = 2 
sr_beta = 2 
np.random.seed(seed=seed_val)
sr_rv = beta.rvs(sr_alpha, sr_beta, size=10000)
np.mean(sr_rv)
sr_p = sr_alpha, sr_beta # sr parameters
plt.figure()
plt.hist(sr_rv, 100, density=True, align='mid')
plt.xlabel('Share of investments')
plt.ylabel('Frequency')
plt.title('Share of investments for reliability ($s_r$), Beta distribution')
plt.text(0.1, 1.7, "$\\alpha$  = "+ str(sr_p[0]))
plt.text(0.1, 1.5, "$\\beta$ = " +str(sr_p[1]))


#Long run demand elasticity, Row 6
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]))

#Displacement of demand
D_alpha = 2 # Range between 0 and 0.75 - skewed towards 0
D_beta = 2
np.random.seed(seed=seed_val)
D_rv = beta.rvs(D_alpha, D_beta, size=1000)
np.mean(D_rv)
D_p = D_alpha, D_beta # D parameters
plt.figure()
plt.hist(D_rv, 100, density=True, align='mid')
plt.xlabel('Share of demand')
plt.ylabel('Frequency')
plt.title('Displacement of demand ($D$), Beta distribution')
plt.text(0.8, 2, "$\\alpha$  = "+ str(D_p[0]))
plt.text(0.8, 1.75, "$\\beta$ = " +str(D_p[1]))

# Lagged effect of displacement on transmission costs
L_alpha = 2
L_beta = 60
np.random.seed(seed=seed_val)
L_rv = beta.rvs(L_alpha,L_beta, size=1000)
np.mean(L_rv)
L_p = L_alpha, L_beta #L parameters
plt.figure()
plt.hist(L_rv, 100, density=True, align='mid')
plt.xlabel('Proportion of cost')
plt.ylabel('Frequency')
plt.title('Lagged effect of displaced demand on costs ($L$), Beta distribution')
plt.text(0.06, 20, "$\\alpha$  = "+ str(L_p[0]))
plt.text(0.06, 17, "$\\beta$ = " +str(L_p[1]))

#Change in generation investment, Row 17
dMW_min = -0.03  
dMW_max = 0 # Range between some max and some min?
np.random.seed(seed=seed_val)
dMW_rv = np.random.uniform(dMW_min,dMW_max, 10000)
dMW_p = dMW_min, dMW_max # dW parameters
plt.figure()
plt.hist(dMW_rv, 100, density=True, align='mid')
plt.xlabel('Proportionate change')
plt.ylabel('Frequency')
plt.title('Change in generation investment ($\\Delta MW$), Uniform distribution')


# MW of generation capacity in export constrained regions (needed for exploring 
# the implications of uncertainty around future baseline investment i.e.
#the basis upon which the previous paremeter is applied).
MW_mu = 4500
MW_sd = 500
np.random.seed(seed=seed_val)
MW_rv = np.random.normal(MW_mu, MW_sd, 10000)
plt.figure()
plt.hist(MW_rv, 100, density=True, align='mid')
MW_p = MW_mu, MW_sd #MW parameters
plt.xlabel('MW in export constrained areas')
plt.ylabel('Frequency')
plt.title('Generation in constrained regions  ($MW$), Normal distribution')
plt.text(2500, 0.0006, "$\\mu$  = "+ str(MW_p[0]))
plt.text(2500, 0.0005, "$\\sigma$ = " +str(MW_p[1]))

    
# Cost benefit function
def f_benefit(QA_p = QA_p, sT_p = sT_p, sr_p = sr_p,
                       dMW_p = dMW_p, eta_p = eta_p, D_p = D_p,
                       MW_p = MW_p, L_p = L_p, ct = ct, Qt = Qt):
    #Parameters
    QA = beta.rvs(QA_p[0], QA_p[1])
    sT = np.random.normal(sT_p[0], sT_p[1])
    sr = beta.rvs(sr_p[0], sr_p[1])
    eta = -1*beta.rvs(eta_p[0], eta_p[1])
    D = beta.rvs(D_p[0], D_p[1])
    L = beta.rvs(L_p[0],L_p[1])
    dMW = np.random.uniform(dMW_p[0],dMW_p[1])
    MW = np.random.normal(MW_p[0], MW_p[1])   
    # Calculations Db = demand benefits, Dc = demand costs, 
    # Genb = generation benefits, Dn = demand net benefits,
    # Nb = Net benefits, Gb = gross benefits
    x = (-1*(eta*((Qt/(QA*Qt)-1)*sT)*Qt*ct*sr*dt))/T #Eq 4
    y = (-1*(eta*((Qt/(QA*Qt)-1)*sT)*Qt*D*L*ct*sr*dt))/T # Eq 6
    z = (-1*dMW*MW*ct*(1-sr)*dt)/T # Eq 7
    Db = np.sum(x)
    Dc = np.sum(y)
    Genb = np.sum(z)
    Dn = Db-Dc
    Nb = Dn+Genb
    Gb = Db+Genb
    
    return Db,Dc,Genb,Dn,Gb,Nb

# Simulation - random
sims = 50000
# Initialise
Db, Dc, Genb, Dn, Gb, Nb = np.zeros(sims), np.zeros(sims), np.zeros(sims), np.zeros(sims), np.zeros(sims), np.zeros(sims)
#Evaluate
for i in range(sims):
    Db[i], Dc[i], Genb[i], Dn[i], Gb[i], Nb[i] = f_benefit()
    
matplotlib.pyplot.close("all")

Nb_trunc = Nb[Nb<np.percentile(Nb,99)]

plt.figure()
plt.hist(Nb_trunc/1000000,100, density=True, align='mid')    
plt.xlabel('Avoided costs ($ millions present value)')
plt.ylabel('Frequency')
plt.title('More efficient generation and demand decisions, monte carlo results')
plt.text(200, 0.0125, 'Mean = %s'%(round(np.mean(Nb)/1000000,2)))

ptiles=[[np.percentile(Nb,5),np.percentile(Nb,50),np.mean(Nb),np.percentile(Nb,95)],
[np.percentile(Gb,5),np.percentile(Gb,50),np.mean(Gb),np.percentile(Gb,95)],
[np.percentile(Genb,5),np.percentile(Genb,50),np.mean(Genb),np.percentile(Genb,95)],
[np.percentile(Db,5),np.percentile(Db,50),np.mean(Db),np.percentile(Db,95)],
[np.percentile(Dc,5),np.percentile(Dc,50),np.mean(Dc),np.percentile(Dc,95)],
[np.percentile(Dn,5),np.percentile(Dn,50),np.mean(Dn),np.percentile(Dn,95)]]

ptiles=pd.DataFrame(ptiles)

ptiles.columns = ["5th percentile","50th percentile","Mean","95th percentile"]
ptiles['Value'] = ['Net benefit', 'Gross benefit','Generation benefit', 'Demand gross benefit', 'Demand cost', 'Demand net benefit']
ptiles.to_csv('Results_demand_and_gen_monte_carlo_final.csv')

Db_trunc = Db[Db<np.percentile(Db,99)]

plt.figure()
plt.hist(Db_trunc/1000000,100, density=True, align='mid')    
plt.xlabel('Avoided costs ($ millions present value)')
plt.ylabel('Frequency')
plt.title('More efficient generation and demand decisions \n demand benefits, monte carlo results')
plt.text(200, 0.0125, 'Mean = %s'%(round(np.mean(Db)/1000000,2)))

Dc_trunc = Dc[Dc<np.percentile(Dc,99)]
plt.figure()
plt.hist(Dc_trunc/1000000,100, density=True, align='mid')    
plt.xlabel('Increased costs ($ millions present value)')
plt.ylabel('Frequency')
plt.title('More efficient generation and demand decisions \n demand costs, monte carlo results')
plt.text(4, 1.25, 'Mean = %s'%(round(np.mean(Dc)/1000000,2)))

plt.figure()
plt.hist(Genb/1000000,100, density=True, align='mid')    
plt.xlabel('Avoided costs ($ millions present value)')
plt.ylabel('Frequency')
plt.title('More efficient generation and demand decisions \n generation benefits, monte carlo results')
plt.text(20, 0.06, 'Mean = %s'%(round(np.mean(Genb)/1000000,2)))