I compare the empirical CDF of a variable with 3 theoretical CDF. I do this for 150 variables and want to print out the result in a single PDF file with 4 charts per page. I do not use a loop but mapply instead. Ideally, I could use par(mfrow=c(2,2)) but I think this works only for R base objects and not ggplot. I looked at gridExtra package here but did not find how to proceed.
library(evd)
library(MASS)
library(fitdistrplus)
library(actuar)
library(ADGofTest)
library (extRemes)
library (lmom)
library(gridExtra)
library(ggplot2)
var1<-rt(10000, df=1)
var2<-rt(10000, df=1)
var3<-rt(10000, df=1)
var4<-rt(10000, df=1)
df<-data.frame(var1,var2, var3, var4)
colnames(df)<-c("var1", "var2", "var3", "var4")
df<-data.frame(var1,var2, var3, var4)
colnames(df)<-c("var1", "var2", "var3", "var4")
pdf()
par(mfrow=c(2,2))
myFUN<-function(x, Name){
empi<-na.omit(x)
empi<-empi[which(empi>0)]
# Theoretical Pareto random series
par.par<-fitdist(empi, "pareto", start=list(shape = 1, scale = 500))
shape.par<-par.par$estimate[1]
scale.par<-par.par$estimate[2]
x.par<-rpareto(NROW(empi), shape.par,scale.par)
# Theoretical Weibull random series
par.wei<-fitdist(empi, "weibull")
shape.wei<-par.wei$estimate[1]
scale.wei<-par.wei$estimate[2]
x.wei<-rweibull(NROW(empi), shape.wei,scale.wei)
# Theoretical GEV random series
# Fittig EVD using the "extRemes" package (can't get it with fitdist)
par.gev <- fevd(empi,type =("GEV"),method=("Lmoments"))
loc.gev<-par.gev$results[1]
shape.gev<-par.gev$results[3]
scale.gev<-par.gev$results[2]
x.gev<-rgev(NROW(empi), loc=loc.gev, scale=scale.gev, shape=shape.gev)
# Create dataframe for using with ggplot+stat_ecdf
df<-data.frame(cbind(empi,rep("Empirical",times=NROW(empi))))
colnames(df)<-c("X","distr")
dfpar<-data.frame(cbind(x.par,rep("Pareto",times=NROW(x.par))))
colnames(dfpar)<-c("X","distr")
dfwei<-data.frame(cbind(x.wei,rep("Weibull",times=NROW(x.wei))))
colnames(dfwei)<-c("X","distr")
dfgev<-data.frame(cbind(x.gev,rep("GEV",times=NROW(x.gev))))
colnames(dfgev)<-c("X","distr")
df<-rbind(df,dfpar)
df<-rbind(df,dfwei)
df<-rbind(df,dfgev)
df$X<-as.numeric(levels(df$X))[df$X]
g<-ggplot(df, aes(X, colour = distr, linetype = distr)) + stat_ecdf(size=1)+theme_classic() +
scale_x_continuous(trans = 'log10')+scale_y_continuous(trans = 'log10') +
xlab("Daily returns")+ylab("CDFs") + ggtitle(Name) + theme(plot.title = element_text(hjust = 0.5)) +
theme(legend.position = c(0.85, 0.25), legend.text=element_text(size=12), legend.title=element_blank())
print(g)
}
allgraph<-mapply(myFUN, df, names(df), SIMPLIFY = FALSE)
dev.off()
Following @bdemarest 's suggestion, I returned to the function gridExtra::marrangeGrob and found a way to do it:
library(evd)
library(MASS)
library(fitdistrplus)
library(actuar)
library(ADGofTest)
library (extRemes)
library (lmom)
library(gridExtra)
library(ggplot2)
var1<-rt(10000, df=1)
var2<-rt(10000, df=1)
var3<-rt(10000, df=1)
var4<-rt(10000, df=1)
df<-data.frame(var1,var2, var3, var4)
colnames(df)<-c("var1", "var2", "var3", "var4")
df<-data.frame(var1,var2, var3, var4)
colnames(df)<-c("var1", "var2", "var3", "var4")
myFUN<-function(x, Name){
empi<-na.omit(x)
empi<-empi[which(empi>0)]
# Theoretical Pareto random series
par.par<-fitdist(empi, "pareto", start=list(shape = 1, scale = 500))
shape.par<-par.par$estimate[1]
scale.par<-par.par$estimate[2]
x.par<-rpareto(NROW(empi), shape.par,scale.par)
# Theoretical Weibull random series
par.wei<-fitdist(empi, "weibull")
shape.wei<-par.wei$estimate[1]
scale.wei<-par.wei$estimate[2]
x.wei<-rweibull(NROW(empi), shape.wei,scale.wei)
# Theoretical GEV random series
# Fittig EVD using the "extRemes" package (can't get it with fitdist)
par.gev <- fevd(empi,type =("GEV"),method=("Lmoments"))
loc.gev<-par.gev$results[1]
shape.gev<-par.gev$results[3]
scale.gev<-par.gev$results[2]
x.gev<-rgev(NROW(empi), loc=loc.gev, scale=scale.gev, shape=shape.gev)
# Create dataframe for using with ggplot+stat_ecdf
df<-data.frame(cbind(empi,rep("Empirical",times=NROW(empi))))
colnames(df)<-c("X","distr")
dfpar<-data.frame(cbind(x.par,rep("Pareto",times=NROW(x.par))))
colnames(dfpar)<-c("X","distr")
dfwei<-data.frame(cbind(x.wei,rep("Weibull",times=NROW(x.wei))))
colnames(dfwei)<-c("X","distr")
dfgev<-data.frame(cbind(x.gev,rep("GEV",times=NROW(x.gev))))
colnames(dfgev)<-c("X","distr")
df<-rbind(df,dfpar)
df<-rbind(df,dfwei)
df<-rbind(df,dfgev)
df$X<-as.numeric(levels(df$X))[df$X]
ggplot(df, aes(X, colour = distr, linetype = distr)) + stat_ecdf(size=1)+theme_classic() +
scale_x_continuous(trans = 'log10')+scale_y_continuous(trans = 'log10') +
xlab("Daily returns")+ylab("CDFs") + ggtitle(Name) + theme(plot.title = element_text(hjust = 0.5)) +
theme(legend.position = c(0.85, 0.25), legend.text=element_text(size=8), legend.title=element_blank())
}
thecharts<-mapply(myFUN, df, names(df), SIMPLIFY = FALSE)
allthecharts<- marrangeGrob(thecharts, nrow=2, ncol=2)
ggsave("allthecharts.pdf", allthecharts)
Previously, I mistakenly added the command g<-ggplot(df,...) print(g) inside the mapply function which gave the error message Error in gList(var1 = list(data = list(list(colour = c("#F8766D", "#F8766D", : only 'grobs' allowed in "gList"