Commit ef42c1e8 authored by Thomas Dimpfl's avatar Thomas Dimpfl
Browse files

add data, literature, R-codes

parent 8066e665

R/.RData

0 → 100644
+262 KiB

File added.

No diff preview for this file type.

View file

R/.Rhistory

0 → 100644
+367 −0
Original line number Diff line number Diff line 
path <- "/home/dimpflth/daten/research/Strom/"

countrylist <- list.files(paste0(path,"/data/volatility_data/"))

countrynames <- str_remove(countrylist, ".csv")

library(stringr)

path <- "/home/dimpflth/daten/research/Strom/"

countrylist <- list.files(paste0(path,"/data/volatility_data/"))

countrynames <- str_remove(countrylist, ".csv")

cn <- 1

indat <- read.csv(paste0(path,"/data/volatility_data/",countrylist[cn]), header=TRUE)

colors <- rep("black", dim(indat)[1])

colors[indat$neg_p>0] <- "red"

y <- as.numeric(substr(indat$date,1,4))

yearlab <- unique(y)

ticks <- c(1,which(diff(y)==1)+1,dim(indat)[1])

labat <- (ticks[2:length(ticks)] + ticks[1:(length(ticks)-1)])/2

plot(sqrt(indat$variance), pch = 20, col=colors, axes=FALSE, xlab="", ylab="volatility")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

plot(indat$daily_span, pch = 20, col=colors, axes=FALSE, xlab="", ylab="Price range")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

summary(sqrt(indat[indat$neg_p==0,"variance"]))

pdf(paste0(path,"plots/descriptive_",countrynames[cn]))

par(mfrow=c(1,2), mar=c(4,4,4,1))

plot(sqrt(indat$variance), pch = 20, col=colors, axes=FALSE, xlab="", ylab="volatility")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

plot(indat$daily_span, pch = 20, col=colors, axes=FALSE, xlab="", ylab="Price range")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

dev.off()

pdf(paste0(path,"plots/descriptive_",countrynames[cn]), width=7, height=4.5)

par(mfrow=c(1,2), mar=c(4,4,4,1))

plot(sqrt(indat$variance), pch = 20, col=colors, axes=FALSE, xlab="", ylab="volatility")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

plot(indat$daily_span, pch = 20, col=colors, axes=FALSE, xlab="", ylab="Price range")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

dev.off()

pdf(paste0(path,"plots/descriptive_",countrynames[cn]), width=7, height=4.5)

par(mfrow=c(1,2), mar=c(4,4,4,1))

plot(sqrt(indat$variance), pch = 20, col=colors, axes=FALSE, xlab="", ylab="volatility")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

plot(indat$daily_span, pch = 20, col=colors, axes=FALSE, xlab="", ylab="Price range")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

mtext(countrynames[cn], outer = TRUE, cex = 1.5, font=4, col=rgb(0.1,0.3,0.5,0.5) )

dev.off()

pdf(paste0(path,"plots/descriptive_",countrynames[cn]), width=7, height=4.5)

par(mfrow=c(1,2), mar=c(4,4,4,1), oma = c(0, 0, 2, 0))

plot(sqrt(indat$variance), pch = 20, col=colors, axes=FALSE, xlab="", ylab="volatility")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

plot(indat$daily_span, pch = 20, col=colors, axes=FALSE, xlab="", ylab="Price range")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

mtext(countrynames[cn], outer = TRUE, cex = 1.5, font=4, col=rgb(0.1,0.3,0.5,0.5) )

dev.off()

pdf(paste0(path,"plots/descriptive_",countrynames[cn]), width=7, height=4.5)

par(mfrow=c(1,2), mar=c(4,4,1,1), oma = c(0, 0, 2, 0))

plot(sqrt(indat$variance), pch = 20, col=colors, axes=FALSE, xlab="", ylab="volatility")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

plot(indat$daily_span, pch = 20, col=colors, axes=FALSE, xlab="", ylab="Price range")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

mtext(countrynames[cn], outer = TRUE, cex = 1.5, font=4, col=rgb(0.1,0.3,0.5,0.5) )

dev.off()

pdf(paste0(path,"plots/descriptive_",countrynames[cn]), width=7, height=4.5)

par(mfrow=c(1,2), mar=c(3,4,1,1), oma = c(0, 0, 2, 0))

plot(sqrt(indat$variance), pch = 20, col=colors, axes=FALSE, xlab="", ylab="volatility")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

plot(indat$daily_span, pch = 20, col=colors, axes=FALSE, xlab="", ylab="Price range")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

mtext(countrynames[cn], outer = TRUE, cex = 1.5, font=4, col=rgb(0.1,0.3,0.5,0.5) )

dev.off()

uhblue <- rgb(0,63,117)

uhblue <- rgb(0,63/255,117/255)

colors <- rep("black", dim(indat)[1])

colors[indat$neg_p>0] <- uhblue

y <- as.numeric(substr(indat$date,1,4))

yearlab <- unique(y)

ticks <- c(1,which(diff(y)==1)+1,dim(indat)[1])

labat <- (ticks[2:length(ticks)] + ticks[1:(length(ticks)-1)])/2

pdf(paste0(path,"plots/descriptive_",countrynames[cn]), width=7, height=4.5)

par(mfrow=c(1,2), mar=c(3,4,1,1), oma = c(0, 0, 2, 0))

plot(sqrt(indat$variance), pch = 20, col=colors, axes=FALSE, xlab="", ylab="volatility")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

plot(indat$daily_span, pch = 20, col=colors, axes=FALSE, xlab="", ylab="Price range")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

mtext(countrynames[cn], outer = TRUE, cex = 1.5, font=4, col=uhblue )

dev.off()

uhblue <- rgb(0,63/255,117/255)

uhred <- rgb(115/225, 13/255, 64/255)

path <- "/home/dimpflth/daten/research/Strom/"

countrylist <- list.files(paste0(path,"/data/volatility_data/"))

countrynames <- str_remove(countrylist, ".csv")

cn <- 1

indat <- read.csv(paste0(path,"/data/volatility_data/",countrylist[cn]), header=TRUE)

colors <- rep("black", dim(indat)[1])

colors[indat$neg_p>0] <- uhred

y <- as.numeric(substr(indat$date,1,4))

yearlab <- unique(y)

ticks <- c(1,which(diff(y)==1)+1,dim(indat)[1])

labat <- (ticks[2:length(ticks)] + ticks[1:(length(ticks)-1)])/2

pdf(paste0(path,"plots/descriptive_",countrynames[cn]), width=7, height=4.5)

par(mfrow=c(1,2), mar=c(3,4,1,1), oma = c(0, 0, 2, 0))

plot(sqrt(indat$variance), pch = 20, col=colors, axes=FALSE, xlab="", ylab="volatility")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

plot(indat$daily_span, pch = 20, col=colors, axes=FALSE, xlab="", ylab="Price range")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

mtext(countrynames[cn], outer = TRUE, cex = 1.5, font=4, col=uhblue )

dev.off()

uhbluebright <- rgb(1,1,0)

colors <- rep("black", dim(indat)[1])

colors[indat$neg_p>0] <- uhbluebright

y <- as.numeric(substr(indat$date,1,4))

yearlab <- unique(y)

ticks <- c(1,which(diff(y)==1)+1,dim(indat)[1])

labat <- (ticks[2:length(ticks)] + ticks[1:(length(ticks)-1)])/2

pdf(paste0(path,"plots/descriptive_",countrynames[cn]), width=7, height=4.5)

par(mfrow=c(1,2), mar=c(3,4,1,1), oma = c(0, 0, 2, 0))

plot(sqrt(indat$variance), pch = 20, col=colors, axes=FALSE, xlab="", ylab="volatility")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

plot(indat$daily_span, pch = 20, col=colors, axes=FALSE, xlab="", ylab="Price range")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

mtext(countrynames[cn], outer = TRUE, cex = 1.5, font=4, col=uhblue )

dev.off()

uhyellow <- rgb(1,204/255,0)

colors <- rep("black", dim(indat)[1])

colors[indat$neg_p>0] <- uhyellow

y <- as.numeric(substr(indat$date,1,4))

yearlab <- unique(y)

ticks <- c(1,which(diff(y)==1)+1,dim(indat)[1])

labat <- (ticks[2:length(ticks)] + ticks[1:(length(ticks)-1)])/2

pdf(paste0(path,"plots/descriptive_",countrynames[cn]), width=7, height=4.5)

par(mfrow=c(1,2), mar=c(3,4,1,1), oma = c(0, 0, 2, 0))

plot(sqrt(indat$variance), pch = 20, col=colors, axes=FALSE, xlab="", ylab="volatility")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

plot(indat$daily_span, pch = 20, col=colors, axes=FALSE, xlab="", ylab="Price range")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

mtext(countrynames[cn], outer = TRUE, cex = 1.5, font=4, col=uhblue )

dev.off()

colors <- rep(uhblue, dim(indat)[1])

colors[indat$neg_p>0] <- uhyellow

y <- as.numeric(substr(indat$date,1,4))

yearlab <- unique(y)

ticks <- c(1,which(diff(y)==1)+1,dim(indat)[1])

labat <- (ticks[2:length(ticks)] + ticks[1:(length(ticks)-1)])/2

pdf(paste0(path,"plots/descriptive_",countrynames[cn]), width=7, height=4.5)

par(mfrow=c(1,2), mar=c(3,4,1,1), oma = c(0, 0, 2, 0))

plot(sqrt(indat$variance), pch = 20, col=colors, axes=FALSE, xlab="", ylab="volatility")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

plot(indat$daily_span, pch = 20, col=colors, axes=FALSE, xlab="", ylab="Price range")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

mtext(countrynames[cn], outer = TRUE, cex = 1.5, font=4, col=uhblue )

dev.off()

## Plots where days with negative prices are highlighted

for (cn in 1:length(countrylist)){

indat <- read.csv(paste0(path,"/data/volatility_data/",countrylist[cn]), header=TRUE)

colors <- rep(uhblue, dim(indat)[1])

colors[indat$neg_p>0] <- uhyellow

y <- as.numeric(substr(indat$date,1,4))

yearlab <- unique(y)

ticks <- c(1,which(diff(y)==1)+1,dim(indat)[1])

labat <- (ticks[2:length(ticks)] + ticks[1:(length(ticks)-1)])/2

pdf(paste0(path,"plots/descriptive_",countrynames[cn]), width=7, height=4)

par(mfrow=c(1,2), mar=c(3,4,1,1), oma = c(0, 0, 2, 0))

plot(sqrt(indat$variance), pch = 20, col=colors, axes=FALSE, xlab="", ylab="volatility")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

plot(indat$daily_span, pch = 20, col=colors, axes=FALSE, xlab="", ylab="Price range")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

mtext(countrynames[cn], outer = TRUE, cex = 1.5, font=4, col=uhblue )

dev.off()

}

cn <- 1

head(indat)

summary(indat[indat$neg_p==0,"variance"])

mean(indat[indat$neg_p==0,"variance"])

mean(indat[indat$neg_p==0,"variance"])

summary((indat[indat$neg_p==0,"daily_span"]))

mean(indat[indat$neg_p==0,"variance"])

summary((indat[indat$neg_p==0,"daily_span"]))

mean(indat[indat$neg_p==0,"variance"])

mean(indat[indat$neg_p==0,"variance"])

mean(indat[indat$neg_p>0,"variance"])

head(indat)

output <- append(output, c(

mean(indat[indat$neg_p==0,"variance"]),

mean(indat[indat$neg_p>0,"variance"]),

mean(indat[indat$neg_p==0,"volatility"]),

mean(indat[indat$neg_p>0,"volatility"]),

mean(indat[indat$neg_p==0,"daily_span"]),

mean(indat[indat$neg_p>0,"daily_span"])

)

output

output <- NULL

output <- append(output, c(

mean(indat[indat$neg_p==0,"variance"]),

mean(indat[indat$neg_p>0,"variance"]),

mean(indat[indat$neg_p==0,"volatility"]),

mean(indat[indat$neg_p>0,"volatility"]),

mean(indat[indat$neg_p==0,"daily_span"]),

mean(indat[indat$neg_p>0,"daily_span"])

)

output

output <- NULL

output <- append(output, c(

mean(indat[indat$neg_p==0,"variance"]),

mean(indat[indat$neg_p>0,"variance"]),

mean(indat[indat$neg_p==0,"volatility"]),

mean(indat[indat$neg_p>0,"volatility"]),

mean(indat[indat$neg_p==0,"daily_span"]),

mean(indat[indat$neg_p>0,"daily_span"])

))

output

## Summary statistics

output <- NULL

for (cn in 1:length(countrylist)){

indat <- read.csv(paste0(path,"/data/volatility_data/",countrylist[cn]), header=TRUE)

output <- append(output, c(

mean(indat[indat$neg_p==0,"variance"]),

mean(indat[indat$neg_p>0,"variance"]),

mean(indat[indat$neg_p==0,"volatility"]),

mean(indat[indat$neg_p>0,"volatility"]),

mean(indat[indat$neg_p==0,"daily_span"]),

mean(indat[indat$neg_p>0,"daily_span"])

))

}

rownames(output) <- countrynames

output

output <- matrix(NA, length(countrylist), 6)

for (cn in 1:length(countrylist)){

indat <- read.csv(paste0(path,"/data/volatility_data/",countrylist[cn]), header=TRUE)

output[cn,] <- c(

mean(indat[indat$neg_p==0,"variance"]),

mean(indat[indat$neg_p>0,"variance"]),

mean(indat[indat$neg_p==0,"volatility"]),

mean(indat[indat$neg_p>0,"volatility"]),

mean(indat[indat$neg_p==0,"daily_span"]),

mean(indat[indat$neg_p>0,"daily_span"])

)

}

rownames(output) <- countrynames

colnames(output) <- c("variance", "variance_n", "volatility", "volatility_n", "daily_span", "daily_span_n")

write.csv(output, "descriptives.csv")

fit <- lm(indat$volatility~indat$neg_p)

fit

## Summary statistics

output <- matrix(NA, length(countrylist), 9)

for (cn in 1:length(countrylist)){

indat <- read.csv(paste0(path,"/data/volatility_data/",countrylist[cn]), header=TRUE)

output[cn,] <- c(

mean(indat[indat$neg_p==0,"variance"]),

mean(indat[indat$neg_p>0,"variance"]),

mean(indat[indat$neg_p>0,"variance"])/mean(indat[indat$neg_p==0,"variance"]),

mean(indat[indat$neg_p==0,"volatility"]),

mean(indat[indat$neg_p>0,"volatility"]),

mean(indat[indat$neg_p>0,"volatility"])/mean(indat[indat$neg_p==0,"volatility"]),

mean(indat[indat$neg_p==0,"daily_span"]),

mean(indat[indat$neg_p>0,"daily_span"]),

mean(indat[indat$neg_p>0,"daily_span"])/ mean(indat[indat$neg_p>0,"daily_span"])

)

}

rownames(output) <- countrynames

colnames(output) <- c("variance", "variance_n", "ratio", "volatility", "volatility_n" "ratio", "daily_span", "daily_span_n", "ratio")

colnames(output) <- c("variance", "variance_n", "ratio", "volatility", "volatility_n", "ratio", "daily_span", "daily_span_n", "ratio")

write.csv(output, "descriptives.csv")

library(stringr)

uhblue <- rgb(0,63/255,117/255)

uhred <- rgb(115/225, 13/255, 64/255)

uhbluebright <- rgb(1,1,0)

uhyellow <- rgb(1,204/255,0)

path <- "/home/dimpflth/daten/research/Strom/"

countrylist <- list.files(paste0(path,"/data/volatility_data/"))

countrynames <- str_remove(countrylist, ".csv")

for (cn in 1:length(countrylist)){

indat <- read.csv(paste0(path,"/data/volatility_data/",countrylist[cn]), header=TRUE)

colors <- rep(uhblue, dim(indat)[1])

colors[indat$neg_p>0] <- uhyellow

y <- as.numeric(substr(indat$date,1,4))

yearlab <- unique(y)

ticks <- c(1,which(diff(y)==1)+1,dim(indat)[1])

labat <- (ticks[2:length(ticks)] + ticks[1:(length(ticks)-1)])/2

pdf(paste0(path,"plots/descriptive_",countrynames[cn]), width=7, height=4)

par(mfrow=c(1,2), mar=c(3,4,1,1), oma = c(0, 0, 2, 0))

plot(sqrt(indat$variance), pch = 20, col=colors, axes=FALSE, xlab="", ylab="volatility")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

plot(indat$daily_span, pch = 20, col=colors, axes=FALSE, xlab="", ylab="Price range")

axis(2)

axis(1, labels=yearlab, at=labat, tick=FALSE)

axis(1, tick=TRUE, at = ticks, labels=FALSE)

box()

mtext(countrynames[cn], outer = TRUE, cex = 1.5, font=4, col=uhblue )

dev.off()

}

R/20240312.R

0 → 100644
+81 −0
Original line number Diff line number Diff line 
###########################################

## Volatility in Power Prices

## Data preparation

###########################################





library(stringr)

library(zoo)



path <- "/home/dimpflth/daten/research/Strom/"



countrylist <- list.files(paste0(path,"/data/hourly_data/"))[-1]

countrynames <- str_remove(countrylist, ".csv")





# compute volatility



for (cn in 1:length(countrylist)){





#     cn <- 1

    indat <- read.csv(paste0(path,"/data/hourly_data/",countrylist[cn]), header=TRUE)

#     indat$date <- as.Date(indat$Datetime..Local.,format="%Y-%m-%d" )

    indat$date <- as.Date(indat$Datetime..UTC.,format="%Y-%m-%d" )

    tmpdate <- unique(indat$date)



    for (td in 1:length(tmpdate)){



        tmpdata <- indat[indat$date == tmpdate[td],]

        n <- dim(tmpdata)[1]

#         tmpvola <- sum( (log(tmpdata[2:n,"Price..EUR.MWhe."]) - log(tmpdata[1:(n-1),"Price..EUR.MWhe."]) )^2)

#         tmpvola <- sum( (  (tmpdata[2:n,"Price..EUR.MWhe."] - tmpdata[1:(n-1),"Price..EUR.MWhe."])/tmpdata[1:(n-1),"Price..EUR.MWhe."] )^2)

#         tmpvola <- sum( (  (tmpdata[2:n,"Price..EUR.MWhe."] - tmpdata[1:(n-1),"Price..EUR.MWhe."]) )^2)

#         tmpvola <- as.numeric(sum(diff(log(as.complex(tmpdata[,"Price..EUR.MWhe."]))))^2)

        tmpvola <- sum((asinh(tmpdata[2:n,"Price..EUR.MWhe."]) - asinh(tmpdata[1:(n-1),"Price..EUR.MWhe."]))^2)   # Schneider Power spot price models with negative prices



        # check if prices were negative

        t_neg_p <- sum(tmpdata[,"Price..EUR.MWhe."]<0)

        t_span <- max(tmpdata[,"Price..EUR.MWhe."]) - min(tmpdata[,"Price..EUR.MWhe."])



        if (td>1){

            output_inner <- rbind(output_inner, c(tmpdate[td],tmpvola,sqrt(tmpvola),t_span,t_neg_p))

        } else {

            output_inner <- data.frame(date=tmpdate[td], variance=tmpvola, volatility=sqrt(tmpvola), daily_span = t_span, neg_p = t_neg_p)

        }

    }



    output_inner$var5 <- c(rep(NA,4), rollapply(output_inner$variance, width=5, sum))

    output_inner$var22 <- c(rep(NA,21), rollapply(output_inner$variance, width=22, sum))

    output_inner$vol5 <- c(rep(NA,4), rollapply(output_inner$volatility, width=5, sum))

    output_inner$vol22 <- c(rep(NA,21), rollapply(output_inner$volatility, width=22, sum))







#     if (cn >1){

#         output <- merge(output, output_inner, by.y = "date", by.x = "date", all.x = TRUE)

#         colnames(output)[cn+1] <- countrynames[cn]

#     } else {

#         output <- output_inner

#         colnames(output)[2] <- countrynames[cn]

#     }



    write.csv(output_inner, paste0(path,"data/volatility_data/",countrylist[cn]))





}



# write.csv(output,"test.csv")

R/20240313.R

0 → 100644
+87 −0
Original line number Diff line number Diff line 
###########################################

## Volatility in Power Prices

## Descriptives

###########################################



library(stringr)



uhblue <- rgb(0,63/255,117/255)

uhred <- rgb(115/225, 13/255, 64/255)

uhbluebright <- rgb(1,1,0)

uhyellow <- rgb(1,204/255,0)



path <- "/home/dimpflth/daten/research/Strom/"



countrylist <- list.files(paste0(path,"/data/volatility_data/"))

countrynames <- str_remove(countrylist, ".csv")





## Plots where days with negative prices are highlighted



for (cn in 1:length(countrylist)){



  indat <- read.csv(paste0(path,"/data/volatility_data/",countrylist[cn]), header=TRUE)

  

  

  colors <- rep(uhblue, dim(indat)[1])

  colors[indat$neg_p>0] <- uhyellow

  

  y <- as.numeric(substr(indat$date,1,4))

  yearlab <- unique(y)

  ticks <- c(1,which(diff(y)==1)+1,dim(indat)[1])

  labat <- (ticks[2:length(ticks)] + ticks[1:(length(ticks)-1)])/2

  

  pdf(paste0(path,"plots/descriptive_",countrynames[cn]), width=7, height=4)

  par(mfrow=c(1,2), mar=c(3,4,1,1), oma = c(0, 0, 2, 0))

  

  plot(sqrt(indat$variance), pch = 20, col=colors, axes=FALSE, xlab="", ylab="volatility")

      axis(2)

      axis(1, labels=yearlab, at=labat, tick=FALSE)

      axis(1, tick=TRUE, at = ticks, labels=FALSE)

      box()

  

  plot(indat$daily_span, pch = 20, col=colors, axes=FALSE, xlab="", ylab="Price range")

      axis(2)

      axis(1, labels=yearlab, at=labat, tick=FALSE)

      axis(1, tick=TRUE, at = ticks, labels=FALSE)

      box()

      

  mtext(countrynames[cn], outer = TRUE, cex = 1.5, font=4, col=uhblue )    

      

  dev.off()    



}





## Summary statistics



output <- matrix(NA, length(countrylist), 9)



for (cn in 1:length(countrylist)){

  

  indat <- read.csv(paste0(path,"/data/volatility_data/",countrylist[cn]), header=TRUE)

  

  output[cn,] <- c(

    mean(indat[indat$neg_p==0,"variance"]),

    mean(indat[indat$neg_p>0,"variance"]),

    mean(indat[indat$neg_p>0,"variance"])/mean(indat[indat$neg_p==0,"variance"]),

    mean(indat[indat$neg_p==0,"volatility"]),

    mean(indat[indat$neg_p>0,"volatility"]),  

    mean(indat[indat$neg_p>0,"volatility"])/mean(indat[indat$neg_p==0,"volatility"]),

    mean(indat[indat$neg_p==0,"daily_span"]),

    mean(indat[indat$neg_p>0,"daily_span"]),

    mean(indat[indat$neg_p>0,"daily_span"])/ mean(indat[indat$neg_p>0,"daily_span"])

    )

}



rownames(output) <- countrynames

colnames(output) <- c("variance", "variance_n", "ratio", "volatility", "volatility_n", "ratio", "daily_span", "daily_span_n", "ratio")

write.csv(output, "descriptives.csv")

R/20240502.R

0 → 100644
+74 −0
Original line number Diff line number Diff line 
###########################################

## Volatility in Power Prices

## Compare transformations

###########################################





library(zoo)



path <- "/home/dimpflth/daten/research/Strom/"



cpdata <- read.csv(paste0(path,"data/comparison_data.csv"), header = TRUE)



colnames(cpdata) <- c("date", "WTI1", "WTI2", "Brent", "CrudeOil", "NaturalGas", "Gasoline", "EPXPowerUK")



# remove negative oil price

cpdata <- cpdata[-which(cpdata$date == "20/04/2020"),]





T <- dim(cpdata)[1]

N <- dim(cpdata)[2]



asret <- asinh(cpdata[2:T,3:N])  - asinh(cpdata[1:(T-1),3:N])

logret <- log(cpdata[2:T,3:N])  - log(cpdata[1:(T-1),3:N])



y <- as.numeric(substr(cpdata$date, 7,10))

yearlab <- unique(y)

ticks <- c(1,which(diff(y)==1)+1,T)

labat <- (ticks[2:length(ticks)] + ticks[1:(length(ticks)-1)])/2



# plot all and compare to normal distribution



pdf(paste0(path,"plots/comparison_data.pdf"))

par(mfrow=c(3,2), mar=c(3,3,1,1))



for (j in 1:(N-2)){







    plot(asret[,j], type="l", col="blue", ylab="", xlab="", axes=FALSE)

    par(new=TRUE)

    plot(logret[,j], type="l", col="gray", axes = FALSE, ylab="", xlab="")

    axis(1, labels=yearlab, at=labat, tick=FALSE)

    axis(1, tick=TRUE, at = ticks, labels=FALSE)

    axis(2)

    box()



    plot(density(asret[,j]), col="blue", main="")

    par(new=TRUE)

    plot(density(logret[,j]), col="gray", main="")





}



dev.off()



stat <- function(func, dat1, dat2){



output <- matrix(NA,length(dat1),3)



    for (j in 1:(length(dat1))){

        output[j,1] <- func(dat1[,j])

        output[j,2] <- func(dat2[,j])

        }

    output[,3] <- output[,2] - output[,1]



        return(output)

}



stat(mean, asret, logret)

stat(median, asret, logret)