Authors: Felipe Van de Sande Araujo, Cristina Pimenta de Mello Spineti Luz, Leonardo Lima Gomes, Luiz Eduardo Teixeira Brandão
Abstract: Despite its continental size and integrated electrical system, Brazil does not have an exchange for trading forward and futures contracts for electricity. Thus, price information for long-term contracts is often obtained through market research and expert opinions. This article proposes a simple yet efficient approach to estimate the forward price of electricity in the Brazilian energy market. The model is based on the equilibrium between two representative agents negotiating bilateral contracts where the agents’ risk aversion is derived from the utility functions related to the Generalized Extended Conditional Value-at-Risk Preference. This model is comprehensive and can be applied to all agents participating in the electricity futures transaction independent of whether they are directly involved in the production chain or simply carry speculative positions. Our results indicate that the model’s forecasted prices, which are based on the participants’ expected behavior, can be used as an indicator for the forward price of electricity, providing more transparency and security for the participants in this market.
This work presents the calculations done for the article referenced above. The following calculations show an optimization of the parameters of an ECP_G risk measure for two agents transacting bilateral contracts in the power market. All the code was run in RStudio using the version of the software below.
R.version _
platform x86_64-w64-mingw32
arch x86_64
os mingw32
system x86_64, mingw32
status
major 4
minor 0.3
year 2020
month 10
day 10
svn rev 79318
language R
version.string R version 4.0.3 (2020-10-10)
nickname Bunny-Wunnies Freak Out # Set the random seed
set.seed(2020) # The year of the start of this projectThe first function is the determination of the ECPG phi which is the certainty equivalent of the agent when it is directly exposed to the spot prices.
ECPG.phi <- function(parameters, dataset, auxList, buyer=TRUE) {
if (any(parameters < 0) || any(parameters > 1)) return(NA)
sumLambda = parameters['Lambda1'] + parameters['Lambda2']
if (sumLambda > 1) return(NA)
Lambda0 = 1 - sumLambda
if (Lambda0 < 0) return(NA)
if (buyer) {
mean.pld = auxList$meanPLD
VaR0 = auxList$VaR0
VaR1 = auxList$VaR1
} else {
dataset = -dataset
mean.pld = -auxList$meanPLD
VaR0 = auxList$VaR1
VaR1 = auxList$VaR0
}
VaRAlpha1 = apply(dataset, 1, quantile, (1-parameters['Alpha1']))
VaRAlpha2 = apply(dataset, 1, quantile, (1-parameters['Alpha2']))
CVaR1 = CVaR2 = NA
CVaR1.data = dataset
CVaR1.data[dataset > VaRAlpha1] <- NA
CVaR2.data = dataset
CVaR2.data[dataset > VaRAlpha2] <- NA
CVaR1 = apply(CVaR1.data, 1, mean, na.rm=TRUE)
CVaR2 = apply(CVaR2.data, 1, mean, na.rm=TRUE)
Sum1 = parameters['Lambda1']*VaRAlpha1 + parameters['Lambda2']*VaRAlpha2
Sum2 = parameters['Lambda1']*(CVaR1 - VaRAlpha1) + parameters['Lambda2']*(CVaR2 - VaRAlpha2)
Limit1 = Lambda0*VaR0 + Sum1
Limit2 = Lambda0*VaRAlpha1 + Sum1
Limit3 = Lambda0*VaRAlpha2 + Sum1
Limit4 = Lambda0*VaR1 + Sum1
ECPG = Lambda0*mean.pld + parameters['Lambda1']*CVaR1 + parameters['Lambda2']*CVaR2
nIter = nrow(dataset)
Lambda1 = rep(parameters['Lambda1'], nIter)
Lambda2 = rep(parameters['Lambda2'], nIter)
if (any(ECPG>=Limit2)) Lambda1[ECPG>=Limit2] = 0
if (any(ECPG>=Limit3)) Lambda2[ECPG>=Limit3] = 0
Q = Lambda0 + Lambda1/(1-parameters['Alpha1']) + Lambda2/(1-parameters['Alpha2'])
A = Sum2 + Lambda0*mean.pld +
Lambda1*VaRAlpha1/(1-parameters['Alpha1']) +
Lambda2*VaRAlpha2/(1-parameters['Alpha2'])
Eq = A/Q
if (any(is.na(Eq))) return(NA)
return(Eq)
}
We will use a list with all the parameters and data which will be called ECPG_Object, in order to facilitate the next function calls. The next function will create the list.
createECPGObject <- function(dataset, discountVector, forwardPrice, buyerParam, sellerParam) {
VaR0 = apply(dataset, 1, quantile, 1)
VaR1 = apply(dataset, 1, quantile, 0)
meanPLD = apply(dataset, 1, mean)
ECPGObject = list(
dataset=dataset,
discountVector=discountVector,
forwardPrice=forwardPrice,
buyerParam=buyerParam,
sellerParam=sellerParam,
VaR0=VaR0,
VaR1=VaR1,
meanPLD=meanPLD
)
return(ECPGObject)
}Then we will define a function to obtain the equilibrium price using the calculated phi for each agent. The distinction between the agents must be made because the signal of the price data changes between them, because of the revenue equations.
ECPG.equilibrium <- function (ECPGobject) {
auxList = list(
VaR0 = ECPGobject$VaR0,
VaR1 = ECPGobject$VaR1,
meanPLD = ECPGobject$meanPLD
)
buyer.phi = ECPG.phi(ECPGobject$buyerParam, ECPGobject$dataset, auxList)
seller.phi = -ECPG.phi(ECPGobject$sellerParam, ECPGobject$dataset, auxList, buyer=FALSE)
equilibrium.prices = (buyer.phi + seller.phi)/2
adjusted.prices = equilibrium.prices/ECPGobject$discountVector
averagePrice = yearly.mean(adjusted.prices)
if (length(averagePrice) != length(ECPGobject$forwardPrice)) stop("Error in ECPG.equilibrium function: Length of fowardPrice vector is different from the length of data divided by 12!")
ECPGobject$buyerPhi = buyer.phi
ECPGobject$sellerPhi = seller.phi
ECPGobject$equilibriumPrices = equilibrium.prices
ECPGobject$adjustedPrices=adjusted.prices
ECPGobject$averagePrice=averagePrice
return(ECPGobject)
}
The next function will be used for the optimization with R command optim.
optim.foo <- function(parameters, ECPGobject) {
param_b = parameters[1:4]
param_s = parameters[5:8]
names(param_b) = names(param_s) = c('Lambda1', 'Lambda2', 'Alpha1', 'Alpha2')
ECPGobject$buyerParam = param_b
ECPGobject$sellerParam = param_s
averagePrice = ECPG.equilibrium(ECPGobject)$averagePrice
return(mean((averagePrice - ECPGobject$forwardPrice)^2))
}The next function is used to calculate the discount vector based on the risk-free rates for each data point.
getDiscountVector <- function(index, modifier, anualRate) {
yearAdj = (11+modifier*12)
monthlyRate = (1+anualRate)^(1/12)-1
completeDiscountVector <- rep(NA, yearAdj)
for (i in 1:yearAdj){
completeDiscountVector[i] <- (1+monthlyRate)^(i-1)
}
getDiscountBounds <- function(index, modifier, completeDiscountVector) {
lower.limit = 1+12*modifier-index
upper.limit = 12+12*modifier-index
return(completeDiscountVector[lower.limit:upper.limit])
}
discountVector = as.numeric(sapply(index, getDiscountBounds, modifier=modifier, completeDiscountVector=completeDiscountVector))
return(discountVector)
}We define an acessory function to calculate the yearly mean using the discount vector.
yearly.mean <- function(x) {
len = length(x)
if (len/12 > len %/%12) stop("Error in yearly.mean function: Length of data is not divisible by 12!")
mean.vec=NULL
iter = len/12
for (i in 1:iter) {
mean.vec[i] = mean(x[(1+(i-1)*12):(i*12)])
}
return(mean.vec)
}Set the initial optimization parameters. They must provide the optimization function with a non-NA return, which is a requirement from optim.
buyer.initial.param = c(Lambda1 = 0, Lambda2 = 0, Alpha1 = 0.5, Alpha2 = 0.95) # Initial Parameters
seller.initial.param = c(Lambda1 = 0, Lambda2 = 0, Alpha1 = 0.5, Alpha2 = 0.95) # Initial ParameterSet the annual discount rate (the economy risk-free rate). A sensitivity analysis for this rate is provided in the sensitivity analysis sector (link).
anualRate = 0.05Set the data points, defined by month and year, to be loaded for model training.
data.index = rbind(
c(7, 2019),
c(8, 2019),
c(9, 2019),
c(10, 2019),
c(11, 2019),
c(12, 2019),
c(1, 2020),
c(2, 2020)
)
# Define the number of years ahead that the forward price is applied to.
seriesA = 1Load the forward price for electricity obtained from DCide Energia.
forwardPrice.data = read.csv(file='./Data/ForwardPrices.csv')Loading future spot price (PLD) data.
data.abbrev = paste(month.abb[data.index[, 1]], data.index[,2], sep="_")
price.data = list()
for (n in 1:8) {
price.data[[data.abbrev[n]]] <- t(read.csv(file=paste('./Data/Data', data.abbrev[n], paste0('A', seriesA, '.csv'), sep="_"), col.names=paste(month.abb, data.index[n,2]+seriesA, sep="-")))
}training.discountVector = getDiscountVector(data.index[,1], seriesA, anualRate)
len = length(training.discountVector)
if (len/12 > len %/%12) stop("Error in discountVector allocation: Length of data is not divisible by 12!")
training.forwardPrice = as.numeric(forwardPrice.data[seriesA, 2:9])
if(len/12 != length(training.forwardPrice)) stop("Error in forwardPrice: Length of data is not compatible with discountVector!")
training.data = Reduce(rbind, price.data)
if (nrow(training.data) != len) stop("Error in trainingDataA2: the number of rows is not compatible with the discountVector!")
training.ECPG_object = createECPGObject(dataset = training.data, buyerParam = buyer.initial.param, sellerParam = seller.initial.param, discountVector = training.discountVector, forwardPrice = training.forwardPrice)
Inputting the data object to the equilibrium function will instantaneously provide the equilibrium price given the selected initial parameters.
training.initialResults = ECPG.equilibrium(training.ECPG_object)
print(training.initialResults$adjustedPrices) Jan.2020 Feb.2020 Mar.2020 Apr.2020 May.2020 Jun.2020 Jul.2020 Aug.2020 Sep.2020 Oct.2020
130.31658 113.69551 100.80909 98.40887 93.20563 97.33694 103.22188 101.35799 100.25743 97.07982
Nov.2020 Dec.2020 Jan.2020 Feb.2020 Mar.2020 Apr.2020 May.2020 Jun.2020 Jul.2020 Aug.2020
92.90089 80.19736 172.92066 154.94408 135.66051 119.61474 116.73274 124.16783 124.92827 123.44964
Sep.2020 Oct.2020 Nov.2020 Dec.2020 Jan.2020 Feb.2020 Mar.2020 Apr.2020 May.2020 Jun.2020
119.94780 115.01348 109.66032 102.89738 168.36718 148.11347 135.61932 115.12322 104.33096 112.05194
Jul.2020 Aug.2020 Sep.2020 Oct.2020 Nov.2020 Dec.2020 Jan.2020 Feb.2020 Mar.2020 Apr.2020
119.01736 113.08721 111.84146 105.98660 101.98258 96.38066 154.02275 154.64323 142.27034 126.64928
May.2020 Jun.2020 Jul.2020 Aug.2020 Sep.2020 Oct.2020 Nov.2020 Dec.2020 Jan.2020 Feb.2020
111.28805 114.83834 120.96265 119.45808 113.87291 111.42977 105.26910 97.62996 189.40988 169.72126
Mar.2020 Apr.2020 May.2020 Jun.2020 Jul.2020 Aug.2020 Sep.2020 Oct.2020 Nov.2020 Dec.2020
170.38496 147.99881 131.35556 131.85121 136.38243 137.53981 134.24549 123.21282 117.93440 105.24825
Jan.2020 Feb.2020 Mar.2020 Apr.2020 May.2020 Jun.2020 Jul.2020 Aug.2020 Sep.2020 Oct.2020
219.86712 189.78008 175.68983 164.38756 144.76853 149.16315 151.72033 148.53520 147.37294 141.21168
Nov.2020 Dec.2020 Jan.2021 Feb.2021 Mar.2021 Apr.2021 May.2021 Jun.2021 Jul.2021 Aug.2021
129.41579 117.97894 82.29963 74.83584 71.34283 67.47625 66.05609 73.48056 78.88422 77.43999
Sep.2021 Oct.2021 Nov.2021 Dec.2021 Jan.2021 Feb.2021 Mar.2021 Apr.2021 May.2021 Jun.2021
75.64360 73.43770 70.61789 68.79611 98.54072 88.91851 83.04969 76.25934 73.06132 80.35147
Jul.2021 Aug.2021 Sep.2021 Oct.2021 Nov.2021 Dec.2021
84.68324 84.70222 83.97916 79.65823 76.24243 71.92401 The object carries all the provided parameters, so it is easy to verify after the execution.
Now the optimization will start. It can take up to 10 minutes running in a single core. Parallelism could be also used if need arise.
## Run the optimization - uses proc.time to measure the execution time
ptm = proc.time()
optim.results = optim(c(buyer.initial.param, seller.initial.param), optim.foo, ECPGobject=training.ECPG_object, method="SANN")
timeSpent = proc.time() - ptm
print(timeSpent) user system elapsed
870.14 9.66 882.42 The optimized parameters are extracted from the optimization results.
## Get ECPG parameters from optimization results
training.ECPG_object$buyerParam = optim.results$par[1:4]
training.ECPG_object$sellerParam = optim.results$par[5:8]The optimal parameters are used to calculate the equilibrium prices. The original object is overwritten because there is no data loss, only new information is added.
## Calculate results from optimized parameters
training.ECPG_object = ECPG.equilibrium(training.ECPG_object)
print(training.ECPG_object$adjustedPrices)Jan.2020 Feb.2020 Mar.2020 Apr.2020 May.2020 Jun.2020 Jul.2020 Aug.2020 Sep.2020 Oct.2020 Nov.2020
164.1603 265.1970 216.7372 189.5391 185.1429 197.4354 215.3346 217.5299 229.9580 255.0911 260.1054
Dec.2020 Jan.2020 Feb.2020 Mar.2020 Apr.2020 May.2020 Jun.2020 Jul.2020 Aug.2020 Sep.2020 Oct.2020
243.2866 200.5650 187.7903 169.9147 253.8775 259.8568 271.9785 288.2995 149.2170 149.5756 139.5760
Nov.2020 Dec.2020 Jan.2020 Feb.2020 Mar.2020 Apr.2020 May.2020 Jun.2020 Jul.2020 Aug.2020 Sep.2020
132.7791 127.7754 190.6412 175.4352 164.9706 143.6773 216.9722 225.1108 147.4742 140.3391 140.3123
Oct.2020 Nov.2020 Dec.2020 Jan.2020 Feb.2020 Mar.2020 Apr.2020 May.2020 Jun.2020 Jul.2020 Aug.2020
130.8242 127.3668 261.6383 178.3329 181.3986 169.6619 154.7897 139.5671 143.7310 150.3345 262.1866
Sep.2020 Oct.2020 Nov.2020 Dec.2020 Jan.2020 Feb.2020 Mar.2020 Apr.2020 May.2020 Jun.2020 Jul.2020
266.0690 137.0331 126.6714 118.7016 220.6390 199.6275 202.9962 176.7106 163.2699 160.8085 165.6458
Aug.2020 Sep.2020 Oct.2020 Nov.2020 Dec.2020 Jan.2020 Feb.2020 Mar.2020 Apr.2020 May.2020 Jun.2020
164.3916 159.6192 150.2491 145.5895 268.1365 251.0878 221.8259 209.8095 198.0010 178.6772 181.2455
Jul.2020 Aug.2020 Sep.2020 Oct.2020 Nov.2020 Dec.2020 Jan.2021 Feb.2021 Mar.2021 Apr.2021 May.2021
184.2241 179.1611 180.6341 172.0481 160.1863 144.3433 225.3101 200.1615 183.8775 179.3444 171.4650
Jun.2021 Jul.2021 Aug.2021 Sep.2021 Oct.2021 Nov.2021 Dec.2021 Jan.2021 Feb.2021 Mar.2021 Apr.2021
185.8318 198.3489 213.6070 202.8853 221.1814 220.0302 202.9710 278.8835 110.5687 257.5283 195.2552
May.2021 Jun.2021 Jul.2021 Aug.2021 Sep.2021 Oct.2021 Nov.2021 Dec.2021
191.8043 206.1048 225.9838 199.7817 236.0811 223.8446 241.5765 238.2001 Finally the results can be printed onscreen.
training.ECPG_object$buyerParam Lambda1 Lambda2 Alpha1 Alpha2
0.25159927 0.09390924 0.32334330 0.69293115 training.ECPG_object$sellerParam Lambda1 Lambda2 Alpha1 Alpha2
0.7466865 0.1712721 0.9912755 0.7641934