DTU31761A2: ReEnergy Participation in Market

library(forecast)
library(gplots)

1. Requirement

Assignment 2 concentrates on participation in electricity markets with a portfolio of renewable energy generation. This portfolio actually consists of a single wind farm located in Western Denmark. Assign- ment 2 requires knowledge of both basic market concepts (day-ahead and balancing stages) and market participation aspects following the lectures, exercise sessions, games and reading material from the first 8 weeks of the course.

2. Data Processing

frameData_1 <- function(){
    dat <- read.table("windpowerforecasts.dat", header = T, sep = ";")
    t_out <- dat[ , 1]
    t_in <- dat[ , 2]
    predStep <- dat[ , 3]
    realMeasure <- dat[ , 4] / 10^3
    predValue <- dat[ , 5] / 10^3
    quan_5 <- dat[ , 6] / 10^3
    quan_10 <- dat[ , 7] / 10^3
    quan_15 <- dat[ , 8] / 10^3
    quan_20 <- dat[ , 9] / 10^3
    quan_25 <- dat[ , 10] / 10^3
    quan_30 <- dat[ , 11] / 10^3
    quan_35 <- dat[ , 12] / 10^3
    quan_40 <- dat[ , 13] / 10^3
    quan_45 <- dat[ , 14] / 10^3
    quan_50 <- dat[ , 15] / 10^3
    quan_55 <- dat[ , 16] / 10^3
    quan_60 <- dat[ , 17] / 10^3
    quan_65 <- dat[ , 18] / 10^3
    quan_70 <- dat[ , 19] / 10^3
    quan_75 <- dat[ , 20] / 10^3
    quan_80 <- dat[ , 21] / 10^3
    quan_85 <- dat[ , 22] / 10^3
    quan_90 <- dat[ , 23] / 10^3
    quan_95 <- dat[ , 24] / 10^3
    quan_100 <- quan_95
    quan_0 <- quan_5
    dat.f <- data.frame(t_out, t_in, predStep, predValue, realMeasure, quan_95, quan_90, quan_85, quan_80, quan_75,
                        quan_70, quan_65, quan_60, quan_55, quan_50, quan_45, quan_40, quan_35, quan_30, quan_25,
                        quan_20, quan_15, quan_10, quan_5, quan_100, quan_0)
    # Clean the data
    t_in.time <- as.numeric(substr(dat.f$t_in, 9, 10))
    dat.f["t_in.time"] <- t_in.time
    dat.f <- dat.f[dat.f$t_in.time == 11, ]
    dat.f <- dat.f[is.element(dat.f$predStep, seq(13, 36)), ]
    dat.f <- dat.f[1:8760, ]
    dat.f$t_in.time <- NULL
    dat.f$t_in <- NULL
    # Label the data
    dat.f["time"] <- rep(seq(0, 23), 365)
    vec_month <- as.numeric(substr(dat.f$t_out, 5, 6))
    dat.f["month"] <- vec_month
    vec_day <- as.numeric(substr(dat.f$t_out, 7, 8))
    dat.f["day"] <- vec_day
    # Return clean data
    dat.f$t_out <- NULL
    return(dat.f)
}
dat.f_q_17 <- frameData_1()
rm(frameData_1)

frameData_2 <- function(str_dat){
    dat <- read.csv(str_dat, header = T, sep = ";")
    seq <- dat[ , 1]
    month <- dat[ , 4]
    day <- dat[ , 5]
    week <- dat[ , 6]
    time <- dat[ , 7]
    d <- dat[ , 8]
    d[is.na(d)] <- 0
    up <- dat[ , 9]
    up[is.na(up)] <- 0
    dw <- dat[ , 10]
    dw[is.na(dw)] <- 0
    dat.f <- data.frame(seq, month, week, day, time, d, up, dw)
    return(dat.f)
}
dat.f_p_16 <- frameData_2("price_16.csv")
dat.f_p_16 <- dat.f_p_16[1: 8784, ]  # 366 days
dat.f_p_17 <- frameData_2("price_17.csv")
dat.f_p_17 <- dat.f_p_17[1: 8760, ]  # 365 days
rm(frameData_2)

Calculate Balancing Price in 2016

#' Function to Calculate Balancing Price
cal_price_b <- function(dat, num_days){
    vec_price_b <- rep(0, length(dat$seq))
    for (i in 1: num_days){
        if (!is.na(dat$d[i])){
            if (abs(dat$up[i] - dat$d[i]) <= abs(dat$dw[i] - dat$d[i])){  
                # dat$d[i] is more next to dat$dw[i], choose the other one
                vec_price_b[i] <- dat$dw[i]
            }
            else if (abs(dat$up[i] - dat$d[i]) >= abs(dat$dw[i] - dat$d[i])){
                vec_price_b[i] <- dat$up[i]
            }
        }
    }
    return(vec_price_b)
}
dat.f_p_16["b"] <- cal_price_b(dat.f_p_16, 8784)
dat.f_p_17["b"] <- cal_price_b(dat.f_p_17, 8760)

3. About Day-Ahead Price

3.1 Box Plot of Day-Ahead Price

sum(dat.f_p_16$d) / 8784

[1] 198.5279

par(bty="n")
plot(factor(dat.f_p_16$month), dat.f_p_16$d, main = "Box Plot of Day-Ahead Price By Months in the Year 2016", 
     xlab = "Month", ylab = "Price")
lines(c(1, 12), c(ave_p_16_d, ave_p_16_d), lty = 2, col = "red", lwd = 3)

legend("bottomleft", inset = .02, legend = "Yearly Average", col = "red", lty = 2, lwd = 3)

par(bty="n")
plot(factor(dat.f_p_16$time), dat.f_p_16$d, main = "Box Plot of Day-Ahead Price By Hours in the Year 2016", 
     xlab = "Hour", ylab = "Price")
lines(c(1, 24), c(ave_p_16_d, ave_p_16_d), lty = 2, col = "red", lwd = 3)

legend("bottomleft", inset = .02, legend = "Yearly Average", col = "red", lty = 2, lwd = 3)

plotBoxPlotTwo <- function(dat.f, ave){
    # Divide the data to clusters 
    vec_index <- rep(0, 8784)
    vec_index[(is.element(dat.f$time, c(23, seq(0, 6)))) & (is.element(dat.f$month, c(10, 11, 12, 1)))] <- 1
    vec_index[(is.element(dat.f$time, seq(7, 18))) & (is.element(dat.f$month, c(10, 11, 12, 1)))] <- 2
    vec_index[(is.element(dat.f$time, seq(19, 22))) & (is.element(dat.f$month, c(10, 11, 12, 1)))] <- 3
    vec_index[(is.element(dat.f$time, c(23, seq(0, 6)))) & (is.element(dat.f$month, seq(2, 5)))] <- 4
    vec_index[(is.element(dat.f$time, seq(7, 18))) & (is.element(dat.f$month, seq(2, 5)))] <- 5
    vec_index[(is.element(dat.f$time, seq(19, 22))) & (is.element(dat.f$month, seq(2, 5)))] <- 6
    dat.f_new <- data.frame(d = dat.f$d, index = vec_index)
    # 
    par(bty = "n")
    plot(factor(dat.f_new$index), dat.f_new$d, main = "Box Plot of Day-Ahead Price By Clusters in the Year 2016", 
         xlab = "Cluster", ylab = "Price", xaxt = "n")
    axis(1, at = seq(1, 7), label = c("others", "Oct-Jan/0-6", "Oct-Jan/7-18", "Oct-Jan/19-22", 
        "Feb-May/0-6", "Feb-May/7-18", "Feb-May/19-22"), las = 0, outer = FALSE)
    lines(c(1, 7), c(ave, ave), lty = 2, col = "red", lwd = 3)
    legend("bottomleft", inset = .02, legend = "Yearly Average", col = "red", lty = 2, lwd = 3)
}
plotBoxPlotTwo(dat.f_p_16, ave_p_16_d)

3.1 Correlation Matrix of Prices over One Day

mat_p_16.h <- matrix(nrow = 366, ncol = 24)
for (i in 1: 24){
    mat_p_16.h[ , i] <- dat.f_p_16[(dat.f_p_16$time == i-1), ]$d
}
# Using function heatmap.2 in library gplots to visualize the correlation matrix
heatmap.2(cor(mat_p_16.h[ , 1: 24]), dendrogram = "none", trace="none", Rowv = F, Colv = F)

3.2 Price Prediction using Expected Value

vec_p_mean.day_16 <- rep(0, 24)
for (i in 0: 23){
    vec_p_mean.day_16[i+1]<- mean(dat.f_p_16[(dat.f_p_16$time == i), ]$d[!is.na(dat.f_p_16[(dat.f_p_16$time == i), ]$d)])
}; rm(i)
plot(seq(0, 23), vec_p_mean.day_16, type = "l", col = "blue", lwd = 2, lty = 1, xaxt = "n", bty = "n",
     main = paste("Expected Day-Ahead Price of All Hours in a Future Day"), xlab = "Hour in one Day", ylab = "Price")
axis(1, at = seq(0, 23), label = seq(1, 24), las = 0, outer = FALSE)

3.3 Predict the Day-Ahead Price using Predicted Wind Power Generation

Update the regression relationship between predicted wind power generation and day-ahead price.

Because there will be more wind power generation in other companies as well, so the price may decrease. The aggregation of renewable generation can impact a lot on the day-ahead clearing price and the effect may be stable.

4. Prediction of Balancing Price

4.1 Regression Relationship between Day-ahead Price and Balancing Price

Predict the Balancing Price using Regression between Predicted Day-Ahead Price

plotPriceRegression <- function(dat.f_p){
    reg_p_16 <- glm(dat.f_p_16$b ~ dat.f_p_16$d)
    # confint(reg_p_16, level = 0.95)
    intConf_p_16 <- predict.lm(reg_p_16, interval = "conf", level = 0.95)
    plot(dat.f_p_16$d, dat.f_p_16$b, col = "blue", yaxt = "n", bty = "n", 
        xlab = "Day-Ahead Price", ylab = "Balancing Price", 
        ylim = c(-500, 2500), main = paste("Regression Relationship Between Day-Ahead Price and Balancing Price"))
    lines(c(-400, 800), reg_p_16$coefficients[1] + reg_p_16$coefficients[2] * c(-400, 800), col = "red", lwd = 2)
    lines(dat.f_p_16$d, intConf_p_16[, 2], lty = 2, col = "black", lwd = 0.5)
    lines(dat.f_p_16$d, intConf_p_16[, 3], lty = 2, col = "black", lwd = 0.5)
    axis(2, at = seq(-500, 2500, by = 500), las = 0, outer = FALSE)
    legend("topleft", inset = .02, legend = c("Observation", "Regression", "95% Confidence Interval"), 
        col = c("blue", "red", "black"), lty = c(NULL, 1, 2), lwd = c(NaN, 2, 0.5), pch = c(1, NaN, NaN))
}
plotPriceRegression(dat.f_p_16)

4.2 Multivariate Regression Model Between Day-Ahead Price and Balancing Price

# Construct matrix for price
mat_p_16.d <- matrix(ncol = 366, nrow = 24)
mat_p_16.b <- matrix(ncol = 366, nrow = 24)
for (i in 1: 24){
    mat_p_16.d[i, ] <- dat.f_p_16[(dat.f_p_16$time == i-1),]$d
    mat_p_16.b[i, ] <- dat.f_p_16[(dat.f_p_16$time == i-1),]$b
}
# Calculate the Regression Coefficient
mat <- matrix(ncol = 25, nrow = 24)
vec_mu <- rep(0, 24)
for (i in 1: 24){
    mat[i, ] <- glm(mat_p_16.b[i, ] ~ mat_p_16.d[1, ] + mat_p_16.d[2, ] + mat_p_16.d[3, ] + 
        mat_p_16.d[4, ] + mat_p_16.d[5, ] + mat_p_16.d[6, ] +
        mat_p_16.d[7, ] + mat_p_16.d[8, ] + mat_p_16.d[9, ] +
        mat_p_16.d[10, ] + mat_p_16.d[11, ] + mat_p_16.d[12, ] +
        mat_p_16.d[13, ] + mat_p_16.d[14, ] + mat_p_16.d[15, ] +
        mat_p_16.d[16, ] + mat_p_16.d[17, ] + mat_p_16.d[18, ] +
        mat_p_16.d[19, ] + mat_p_16.d[20, ] + mat_p_16.d[21, ] +
        mat_p_16.d[22, ] + mat_p_16.d[23, ] + mat_p_16.d[24, ])$coefficients[1:25]
}
vec_mu <- mat[ , 1]
mat_c <- mat[ , 2: 25]
rm(mat)
# Check the result
all((mat_c - cov(t(mat_p_16.b), t(mat_p_16.d)) %*% solve(cov(t(mat_p_16.d)))) < 0.000001)

[1] TRUE

5 Trading Strategies

#' Function to trade in day-ahead market and calculate the revenue
tradeDayAhead <- function(vec_e_d, vec_p_d){
    vec_r_d <- vec_e_d * vec_p_d
    return(vec_r_d)
}
#' Function to get the real measured values in that day
getRealValue <- function(vec_p_d_all, vec_p_b_all, vec_e_all, day){
    # list_realValue <- getRealValue(dat.f_p$d, dat.f_p$b, dat.f_q$realMeasure)
    vec_p_d <- vec_p_d_all[(24 * (day-1) + 1): (24 * (day-1) + 24)]
    vec_p_b <- vec_p_b_all[(24 * (day-1) + 1): (24 * (day-1) + 24)]
    vec_e <- vec_e_all[(24 * (day-1) + 1): (24 * (day-1) + 24)]
    return(list(d = vec_p_d, b = vec_p_b, e = vec_e))
}
#' Function to trade in balancing market and calculate the revenue
tradeBalancing <- function(vec_e_d, vec_e, vec_p_up, vec_p_dw){
    vec_e_up <- rep(0, 24)
    vec_e_dw <- rep(0, 24)
    vec_r_b <- rep(0, 24)
    # If na, we don't trade that day anyway
    vec_e[is.na(vec_e)] <- 0
    for (i in 1: 24){
        # 1. Calculate regulation quantity
        if (vec_e[i] <= vec_e_d[i]){  # Output too less
            vec_e_up[i] <- vec_e[i] - vec_e_d[i]
            vec_e_dw[i] <- 0
        } else {  # Output too much
            vec_e_up[i] <- 0
            vec_e_dw[i] <- vec_e[i] - vec_e_d[i]
        }
        # 2. Calculate the revenue in balancing market
        if (vec_e_up[i] > 0){
            print("Error:", vec_e_up[i], "is greater than 0.")
        }
        vec_r_b[i] <- vec_p_up[i] * vec_e_up[i] + vec_p_dw[i] * vec_e_dw[i]
    }
    return(vec_r_b)
}

5.1 Trading Strategy using Deterministic Wind Power Forecasts

5.2 Trading Strategy using Probabilistic Wind Power Forecasts

#' Calculate up-regulation and down-regulation prices using day-ahead price and balancing price
calReguPrice <- function(vec_p_d, vec_p_b){
    # If Na, we don't trade that day anyway
    vec_p_d[is.na(vec_p_d)] <- 0
    vec_p_b[is.na(vec_p_b)] <- 0
    # 
    vec_p_up <- rep(0, 24)
    vec_p_dw <- rep(0, 24)
    for (i in 1: 24){
        # 1. Calculate the predicted regulation price
        if (vec_p_b[i] >= vec_p_d[i]){  # higher balancing price
            vec_p_up[i] <- vec_p_b[i]
            vec_p_dw[i] <- vec_p_d[i]
        } else {  # lower balancing price
            vec_p_up[i] <- vec_p_d[i]
            vec_p_dw[i] <- vec_p_b[i]
        }
    }
    return(list(up = vec_p_up, dw = vec_p_dw))
}
#' Function to calculate the trading volume in day-ahead market using Prob-Fore strategy
useStrategyProb <- function(vec_p_d_hat, vec_p_up_hat, vec_p_dw_hat, dat.f){
    vec_psi_up_hat <- rep(0, 24)
    vec_psi_dw_hat <- rep(0, 24)
    vec_ratio <- rep(0, 24)
    vec_e_d <- rep(0, 24)
    for (i in 1: 24){
        vec_psi_up_hat[i] <- vec_p_up_hat[i] - vec_p_d_hat[i]
        vec_psi_dw_hat[i] <- vec_p_d_hat[i] - vec_p_dw_hat[i]
        vec_ratio[i] <- vec_psi_dw_hat[i] / (vec_psi_up_hat[i] + vec_psi_dw_hat[i])
        str <- paste("quan_", as.integer(vec_ratio[i] * 100 / 5) * 5, sep = "")
        vec_e_d[i] <- dat.f[str][i,1]
    }
    return(vec_e_d)
}
#'
plotPredCDFun <- function(dat.f, day = 1, h = 1){
    seq_prob <- seq(0, 100, by = 5)
    vec_q <- rep(0, length(seq_prob))
    for (i in 1: length(seq_prob)){
        str <- paste("quan_", seq_prob[i], sep = "")
        vec_q[i] <- dat.f[str][(24 * (day - 1) + h),1]
    }
    plot(c(vec_q, 160), c(seq_prob, 100), type = "l", col = "blue", lwd = 2,
         main = paste("CDF and Expected Value of Predicted Wind Power Output in", day, "day at", h, "O'clock"), 
         xlab = "Wind Power Output (MW)", ylab = "Probability", xaxt = "n", yaxt = "n", bty = "n")
    axis(1, at = seq(0, 160, by = 20), las = 0, outer = FALSE)
    axis(2, at = seq(0, 100, by = 20), labels = seq(0, 1, by = 0.2), las = 0, outer = FALSE)
    legend("topleft", inset = .02, legend = c("CDFunction", "Expected Value"), 
            col = c("blue", "red"), lty = c(1, 2), lwd = c(2, 2))
    # Add the expected value in the plot
    predValue <- dat.f$predValue[(24 * (day - 1) + h)]
    lines(c(predValue, predValue), c(0, 100), col = "red", lwd = 2, lty = 2)
    # points(predValue, 50, col = "red", lwd = 2, pch = 4, cex = 3)
}

plotPredCDFun(dat.f_q_17, 114, 19)

6 Performance Evaluation

#' Function to plot the accumulated revenue
#' If vec_rAccuMax.day is not NaN, plot the maximum accumulated revenue as well
plotRevenueDay <- function(vec_rAccu.day, vec_rAccuMax.day = NaN, str_strategy, num_day = 365, year = 2017){
    if (!(is.na(vec_rAccuMax.day[1]))){
        yMin <- min(vec_rAccu.day, vec_rAccuMax.day)
        yMax <- max(vec_rAccu.day, vec_rAccuMax.day)
    } else {
        yMin <- min(vec_rAccu.day)
        yMax <- max(vec_rAccu.day)
    }
    plot(seq(1: num_day), vec_rAccu.day, type = "l", ylim = c(yMin, yMax), col = "blue", lwd = 2,
            main = paste("Accumulated Revenue Over the Year", year, "using", str_strategy, "Strategy"), 
            xlab = "month", ylab = "revenue (DKK)", xaxt = "n", bty = "n")
    # text(365, vec_rAccu.day[365], paste(vec_rAccu.day[365]), pos = 1)
    axis(1, at = c(1, 32, 60, 91, 122, 152, 183, 214, 244, 275, 305, 336),
         substr(month.name, 1, 3), las = 0, outer = FALSE)
    if (!(is.na(vec_rAccuMax.day[1]))){
        legend("topleft", inset = .02, legend = c("Strategic Trading", "Perfect Info"), 
            col = c("blue", "red"), lty = c(1, 2), lwd = c(2, 2))
        lines(seq(1: num_day), vec_rAccuMax.day, type = "l", col = "red", lwd = 2, lty = 2)
    } else {
        legend("topleft", inset = .02, legend = c("Strategy"), 
            col = c("blue"), lty = c(2), lwd = c(2))
    } 
}

6.1 The maximum revenue from trade with Perfect Information:

calRevenueMax <- function(vec_p_d, vec_qReal){
    # If any data is missing, don't trade on that day
    vec_qReal[is.na(vec_qReal)] <- 0
    vec_qReal[is.na(vec_p_d)] <- 0
    vec_p_d[is.na(vec_qReal)] <- 0
    vec_p_d[is.na(vec_p_d)] <- 0
    # Calculate vector of maximum revenue and maximum accumulated revenue
    vec_rMax.h <- vec_p_d * vec_qReal
    vec_rMax.day <- rep(0, 365)
    vec_rAccuMax.day <- rep(0, 365)
    for (day in 1: 365){
        vec_rMax.day[day] <- sum(vec_rMax.h[(24 * (day-1) + 1): (24 * (day-1) + 24)])
        if (day == 1){
            vec_rAccuMax.day[day] <- vec_rMax.day[1]
        } else {
            vec_rAccuMax.day[day] <- vec_rMax.day[day] + vec_rAccuMax.day[day-1]
        }
    }; rm(day)
    plot(seq(1: 365), vec_rAccuMax.day, type = "l", col = "red", lty = 2, lwd = 2,
            main = paste("Maximum Accumulated Revenue Over the Year 2017 with Perfect Information"), 
            xlab = "month", ylab = "revenue (DKK)", xaxt = "n", bty = "n")
    axis(1, at = c(1, 32, 60, 91, 122, 152, 183, 214, 244, 275, 305, 336),
         substr(month.name, 1, 3), las = 0, outer = FALSE)
    return(vec_rAccuMax.day)
}
vec_rAccuMax.day <- calRevenueMax(dat.f_p_17$d, dat.f_q_17$realMeasure)

rm(calRevenueMax)

6.2 Base Trading Strategy using Fixed Volume

vec_rAccu.day <- rep(0, 365)
vec_r.day <- rep(0, 365)
# For every day in 2017, calculate the revenue
for (day in 1: 365){
    # 1. Calculate trade volume
    vec_e_d <- rep(150, 24)
    # If any data is missing, don't trade on that day
    list_real <- getRealValue(dat.f_p_17$d, dat.f_p_17$b, dat.f_q_17$realMeasure, day)
    vec_e_d[is.na(list_real$d)] <- 0
    vec_e_d[is.na(list_real$b)] <- 0
    vec_e_d[is.na(list_real$e)] <- 0
    vec_e_d[(is.na(vec_e_d))] <- 0
    list_real$e[is.na(list_real$d)] <- 0
    list_real$e[is.na(list_real$b)] <- 0
    list_real$e[is.na(list_real$e)] <- 0
    # 2. Trade in day-ahead market
    vec_r_d <- tradeDayAhead(vec_e_d, list_real$d)
    # 3. Trade in balancing marekt
    list_pRegu <- calReguPrice(list_real$d, list_real$b)
    vec_r_b <- tradeBalancing(vec_e_d, list_real$e, list_pRegu$up, list_pRegu$dw)
    # 4. Calculate the revenue in this day
    vec_r.day[day] <- sum(vec_r_d + vec_r_b)
    if (day == 1){
        vec_rAccu.day[day] <- vec_r.day[day]
    } else {
        vec_rAccu.day[day] <- vec_r.day[day] + vec_rAccu.day[day-1]
    }
}; rm(vec_e_d, vec_r_d, list_real, list_pRegu, vec_r_b)
plotRevenueDay(vec_rAccu.day, vec_rAccuMax.day, "\"Fixed Volume\"")

print(paste("Result: total revenue =", vec_rAccu.day[365]))

[1] "Result: total revenue = 129705990.68789"

print(paste("Result: gamma =", vec_rAccu.day[365] / vec_rAccuMax.day[365]))

[1] "Result: gamma = 0.890662217817209"

rm(vec_rAccu.day, vec_r.day)

6.3 Trading Strategy using Deterministic Wind Power Forecasts

vec_rAccu.day <- rep(0, 365)
vec_r.day <- rep(0, 365)
# For every day in 2017, calculate the revenue
for (day in 1: 365){
    # 1. Calculate trade volume
    vec_e_d <- dat.f_q_17$predValue[(24 * (day-1) + 1): (24 * (day-1) + 24)]
    vec_e_d <- 1.0 * vec_e_d
    # If any data is missing, don't trade on that day
    list_real <- getRealValue(dat.f_p_17$d, dat.f_p_17$b, dat.f_q_17$realMeasure, day)
    vec_e_d[is.na(list_real$d)] <- 0
    vec_e_d[is.na(list_real$b)] <- 0
    vec_e_d[is.na(list_real$e)] <- 0
    vec_e_d[(is.na(vec_e_d))] <- 0
    list_real$e[is.na(list_real$d)] <- 0
    list_real$e[is.na(list_real$b)] <- 0
    list_real$e[is.na(list_real$e)] <- 0
    # 2. Trade in day-ahead market
    vec_r_d <- tradeDayAhead(vec_e_d, list_real$d)
    # 3. Trade in balancing marekt
    list_pRegu <- calReguPrice(list_real$d, list_real$b)
    vec_r_b <- tradeBalancing(vec_e_d, list_real$e, list_pRegu$up, list_pRegu$dw)
    # 4. Calculate the revenue in this day
    vec_r.day[day] <- sum(vec_r_d + vec_r_b)
    if (day == 1){
        vec_rAccu.day[day] <- vec_r.day[day]
    } else {
        vec_rAccu.day[day] <- vec_r.day[day] + vec_rAccu.day[day-1]
    }
}; rm(vec_e_d, vec_r_d, list_real, list_pRegu, vec_r_b)
plotRevenueDay(vec_rAccu.day, vec_rAccuMax.day, "\"Deterministic Forecast\"")

print(paste("Result: total revenue =", vec_rAccu.day[365]))

[1] "Result: total revenue = 140634479.24337"

print(paste("Result: gamma =", vec_rAccu.day[365] / vec_rAccuMax.day[365]))

[1] "Result: gamma = 0.965705720454151"

rm(vec_rAccu.day, vec_r.day)

6.4 Trading Strategy using Probablistic Wind Power Forecasts

vec_rAccu.day <- rep(0, 365)
vec_r.day <- rep(0, 365)
# 1. Predict day-ahead and balancing prices
vec_p_d_hat <- vec_p_mean.day_16
vec_p_b_hat <- (mat_c %*% vec_p_d_hat)[ , 1]  # vec_p_b_hat <- rep(0, 24)
list_pRegu_hat <- calReguPrice(vec_p_d_hat, vec_p_b_hat)
# For every day in 2017, calculate the revenue
for (day in 1: 365){
    # 2. Calculate trade volume
    vec_e_d <- useStrategyProb(vec_p_d_hat, list_pRegu_hat$up, list_pRegu_hat$dw, dat.f_q_17)
    # If any data is missing, don't trade on that day
    list_real <- getRealValue(dat.f_p_17$d, dat.f_p_17$b, dat.f_q_17$realMeasure, day)
    vec_e_d[is.na(list_real$d)] <- 0
    vec_e_d[is.na(list_real$b)] <- 0
    vec_e_d[is.na(list_real$e)] <- 0
    list_real$e[is.na(list_real$d)] <- 0 
    list_real$e[is.na(list_real$b)] <- 0
    list_real$e[is.na(list_real$e)] <- 0
    # 3. Trade in day-ahead market
    vec_r_d <- tradeDayAhead(vec_e_d, list_real$d)
    # 4. Trade in balancing marekt
    list_pRegu <- calReguPrice(list_real$d, list_real$b)
    vec_r_b <- tradeBalancing(vec_e_d, list_real$e, list_pRegu$up, list_pRegu$dw)
    # 5. Calculate the revenue in this day
    vec_r.day[day] <- sum(vec_r_d + vec_r_b)
    if (day == 1){
        vec_rAccu.day[day] <- vec_r.day[day]
    } else {
        vec_rAccu.day[day] <- vec_r.day[day] + vec_rAccu.day[day-1]
    }
}; rm(vec_e_d, vec_r_d, list_real, list_pRegu, vec_r_b)
plotRevenueDay(vec_rAccu.day, vec_rAccuMax.day, "\"Prob Forecast\"")

print(paste("Result: total revenue =", vec_rAccu.day[365]))

[1] "Result: total revenue = 133768823.07183"

print(paste("Result: gamma =", vec_rAccu.day[365] / vec_rAccuMax.day[365]))

[1] "Result: gamma = 0.918560785049983"

rm(vec_rAccu.day, vec_r.day)

vec_rAccu.day <- rep(0, 365)
vec_r.day <- rep(0, 365)
# 1. Predict day-ahead and balancing prices
vec_p_d_hat <- vec_p_mean.day_16
vec_p_b_hat <- (mat_c %*% vec_p_d_hat)[ , 1]  # vec_p_b_hat <- rep(0, 24)
list_pRegu_hat <- calReguPrice(vec_p_d_hat, vec_p_b_hat)
# For every day in 2017, calculate the revenue
for (day in 1: 365){
    # 2. Calculate trade volume
    vec_e_d <- useStrategyProb(vec_p_d_hat, list_pRegu_hat$up, list_pRegu_hat$dw, dat.f_q_17)
    vec_e_d <- vec_e_d * 0.7
    # If any data is missing, don't trade on that day
    list_real <- getRealValue(dat.f_p_17$d, dat.f_p_17$b, dat.f_q_17$realMeasure, day)
    vec_e_d[is.na(list_real$d)] <- 0
    vec_e_d[is.na(list_real$b)] <- 0
    vec_e_d[is.na(list_real$e)] <- 0
    list_real$e[is.na(list_real$d)] <- 0 
    list_real$e[is.na(list_real$b)] <- 0
    list_real$e[is.na(list_real$e)] <- 0
    # 3. Trade in day-ahead market
    vec_r_d <- tradeDayAhead(vec_e_d, list_real$d)
    # 4. Trade in balancing marekt
    list_pRegu <- calReguPrice(list_real$d, list_real$b)
    vec_r_b <- tradeBalancing(vec_e_d, list_real$e, list_pRegu$up, list_pRegu$dw)
    # 5. Calculate the revenue in this day
    vec_r.day[day] <- sum(vec_r_d + vec_r_b)
    if (day == 1){
        vec_rAccu.day[day] <- vec_r.day[day]
    } else {
        vec_rAccu.day[day] <- vec_r.day[day] + vec_rAccu.day[day-1]
    }
}; rm(vec_e_d, vec_r_d, list_real, list_pRegu, vec_r_b)
plotRevenueDay(vec_rAccu.day, vec_rAccuMax.day, "\"Prob Forecast\"")

print(paste("Result: total revenue =", vec_rAccu.day[365]))

[1] "Result: total revenue = 133651552.544013"

print(paste("Result: gamma =", vec_rAccu.day[365] / vec_rAccuMax.day[365]))

[1] "Result: gamma = 0.917755514392583"

rm(vec_rAccu.day, vec_r.day)