In paul4forest/GFPMoutput: This program reads output simulation of the Global Forest Product Model (GFPM).
library(knitr)
opts_knit$set(root.dir="../..") # file paths are relative to the root of the project directory
library(GFPMoutput)
library(dplyr)
library(tidyr)
library(ggplot2)
# Attention do not load the tradeflow package here,
# It's clean function masks this one.
selectedproducts <- c("Newsprint", "PWPaper",
"OthPaper", "Sawnwood",
"ParticleB", "FiberB")
Blog
- Before July first try two different base files
the default 2016 world file and base105 from previous
work with Joseph.
- 20160706 try several ways to run the simulation
It seems there is an error after the first period
when the Excel files where modified in a French version of
Excel. I changed the region and Language settings to
English US.
GFPM also gives also an error when there are
positive price elasticities. I decided to only export
those products which have negative price elasticities.
The issue is that simulation results are exactly the same.
Load demand elasticities from the GFPM base scenario
worlddemand <- read.csv("rawdata/World2016demand.csv",
stringsAsFactors = FALSE) %>%
left_join(countrycodes[c("Country_Code","Country")]) %>%
left_join(productcodes)
worlddemand105 <- read.csv("rawdata/Worlddemand105.csv",
stringsAsFactors = FALSE) %>%
left_join(countrycodes[c("Country_Code","Country")]) %>%
left_join(productcodes)
# Country with zero elast
worlddemand %>% filter(priceelast == 0)
message("Dummy country zy can be ignored")
#
elastsummmary <- worlddemand %>%
filter(priceelast != 0) %>%
select(Product, priceelast, gdpelast) %>%
distinct()
elastsummmary %>% kable()
elastsummmary %>%
filter(Product %in% selectedproducts & priceelast !=-0.25) %>%
dput()
Show products that have more than one elasticity
elastsummmary %>% filter(duplicated(Product)) %>% kable()
Show the group of countries which share a same elasticity
of fuelwood demand.
fwdgroup <- worlddemand %>%
filter(Product == "Fuelwood" & !is.na(Country)) %>%
group_by(priceelast,gdpelast) %>%
summarise(Country = paste(Country,collapse=","))
fwdgroup %>% kable()
plygroup <- worlddemand %>%
filter(Product == "Plywood" & !is.na(Country)) %>%
group_by(priceelast, gdpelast) %>%
summarise(Country = paste(Country, collapse=","))
plygroup %>% kable()
newsgroup <- worlddemand %>%
filter(Product == "Newsprint" & !is.na(Country)) %>%
group_by(priceelast, gdpelast) %>%
summarise(Country = paste(Country, collapse=","))
newsgroup %>% kable()
Update the elasticities table for the World.xls of GFPM
Use elasticities from the pannel cointegration (PC) article .
productnames <- data_frame(Product = c("Newsprint", "PWPaper", "OthPaper", "Sawnwood", "ParticleB", "FiberB", "Plywood"),
item = c("Newsprint", "PrintingandWritingPaper", "OtherPaperandPaperboard",
"SawnwoodConiferous", "ParticleBoard", "Fibreboard", "Plywood")) %>%
mutate(item = factor(item, levels = item, ordered = TRUE))
# GDP, Price as in the article's table
elast <- read.csv("/home/paul/R/forestproductsdemand/data-end/elastforGFPM.csv", stringsAsFactors = FALSE) %>%
rename(gdpelastpc = lgdpreur,
priceelastpc = lpricereur) %>%
# Rename Items to GFPM Products
left_join(productnames, by = "item") %>%
select(-item)
elast %>% kable()
# Simangunsong product name missing before sawnwood
# not important Total swd and total panel not useful anyway
# Only plywood would be interesting.
# Compare with elasticites currently in the GFPM
elast %>%
left_join(filter(elastsummmary, Product %in% elast$Product), by="Product") %>%
select(Product, estimator, gdpelast, gdpelastpc, priceelast, priceelastpc) %>%
kable()
# Chunk set to eval=FALSE because it's used only once before a simulation
elast %>% filter(estimator == "dols") %>%
adddemandelast(worlddemand105) %>%
write.csv("rawdata/World105demanddols.csv")
elast %>% filter(estimator == "pmg") %>%
adddemandelast(worlddemand105) %>%
write.csv("rawdata/World105demandpmg.csv")
elast %>% filter(estimator == "dols") %>%
adddemandelast(worlddemand) %>%
write.csv("rawdata/World2016demanddols.csv")
elast %>% filter(estimator == "pmg") %>%
adddemandelast(worlddemand) %>%
write.csv("rawdata/World2016demandpmg.csv")
elast %>% filter(estimator == "simangunsong") %>%
adddemandelast(worlddemand) %>%
write.csv("rawdata/World2016demandsimangunsong.csv")
if(FALSE){
# Only Sawnwood
message("DOLS::Newsprint and Sawnwood have a positive price elasticity")
message("PMG::Newsprint has a positive price elasticity")
# DOLS elast without Newsprint without sawnwood
elast %>%
filter(estimator == "dols" &
! Product %in% c("Newsprint","Sawnwood")) %>%
adddemandelast(worlddemand) %>%
write.csv("rawdata/World2016demandnonewsnoswd.csv")
# Only particle board
elast %>% filter(estimator == "dols" &
Product == "ParticleB") %>%
adddemandelast(worlddemand) %>%
write.csv("rawdata/World2016demanddolsParticleB.csv")
# Only Othpaper
elast %>% filter(estimator == "dols" &
Product == "OthPaper") %>%
adddemandelast(worlddemand) %>%
write.csv("rawdata/World2016demanddolsOthPaper.csv")
# Pmg elast without newsprint,
# Avoid newsprint because it's positive price elasticity
# seems to be creating an error
elast %>% filter(estimator == "pmg" &
Product != "Newsprint") %>%
adddemandelast(worlddemand) %>%
write.csv("rawdata/World2016demandpmgnonews.csv")
}
if (FALSE){
pcelast <- elast %>% filter(estimator == "dols")
# Develop the adddemandelast function
worlddemand2 <- worlddemand %>%
left_join(pcelast, by="Product")
# Replace those that have a new elasticity by the new value
worlddemand2 <- worlddemand2 %>%
filter(!is.na(gdpelastpc)) %>%
mutate(gdpelast = gdpelastpc,
priceelast = priceelastpc) %>%
# Add back those that didn't have a new elasticity
rbind(filter(worlddemand2,is.na(gdpelastpc))) %>%
arrange(Country_Code, Product_Code)
# Did the number of rows stay the same?
stopifnot(identical(nrow(worlddemand),nrow(worlddemand2)))
# Did the order of rows stay the same?
stopifnot(identical(worlddemand[c("Country_Code","Product_Code","basedemand")],worlddemand2[c("Country_Code","Product_Code","basedemand")]))
write.csv(worlddemand2, "rawdata/World2016demandPC.csv",
row.names = FALSE)
worlddemand2016dols <- elast %>%
filter(estimator == "dols") %>%
adddemandelast(worlddemand)
}
GFPM issues
Positive price elasticities
Positive price elasticities cause the GFPM simulation to break after the first optimisation round.
- The simulation breaks when I use all PMG estimates including newsprint elasticitity of demand.
- The simulation runs find when I use PMG estimates excluding newsprint elasticity of demand
- Similar issue for DOLS estimatesof the newsprint and sawnwood elasticity of demand.
Therefore I used only negative price elasticities.
Price elasticities should be negative according to the theory.
Remove Exogchange elasticities changes
Joseph Buongiorno:
The reason for the small difference that you got is that in the world file that you used, the elasticities changed over time. This is reflected in the ExogChange file where the elasticities differ in each period. It looks like you have not changed that. So, the results reflect only the effect of a change in elasticity in the first period, which is small. I suggest changing the elasticities in all periods in ExogChange so that they are the same as in the first period and equal to your revised elasticities.
I removed GDP and price elasticity from exogenous change.
By setting the GDP elasticity (column R) to 0 with a double click.
Double click copies a full column down untill the next
empty cell.
Then I copied this column of 0 to the price elasticity (column H).
Otherwise a double click on the price elastcitity section
would errase all content from column H, for all data types.
WHich is not what I want. I wand to modify the D data types only.
Excel issues
Because I use a French/ English version of excel
I had to change dots to commas in the file.
I swidched windows regional and language settings to US so that
I don't have to do this commas / dot replacement each time I load a csv file.
Then here is how to send updated elasticities to GFPM:
- open this csv in Excel, convert the text to columns if
necessary using comma as a separator.
- paste the priceelast and gdpelast columns in the
World.xls table
- replace points by commas because we are dealing with
a french version of excel
(Side note: commas are incompatible with csv, they would require every field to be quoted. Now I understand why read.csv2 accepts text separated by semi-colons)
Then run the simulation as an alternative scenario
and go back to this document
to load both the base and alternative scenarios below.
Load simulation data
I copied the file to an archive called "pelpsbase2016"
under ./rawdata.
This needs to be evaluated only once.
Once data has been processed and storred in enddata.
It can be read from there.
message("Chunk set to eval=FALSE as it needs to be evaluated only once to load data from GFPM after a simulation.")
savePELPSToRdata("pelpsbase2016", compressed = "zip")
savePELPSToRdata("pelpspc2016", compressed = "zip")
savePELPSToRdata("pelpslow2016", compressed = "zip")
base2016 <- clean("pelpsbase2016.RDATA","base2016")
pc <- clean("pelpspc2016.RDATA","pc2016")
low <- clean("pelpslow2016.RDATA","low2016")
# Write to endata
message("Use bindScenarios() to combine scenarios ",
"in one list, then write the two scenarios in one rds file.")
scenarios <- Reduce(bindScenarios,
list(base, pc, low))
saveRDS(scenarios, "enddata/baseandpc2016.rds")
# base105 scenarios
savePELPSToRdata("base105low", compressed = "zip")
savePELPSToRdata("base105high", compressed = "zip")
base105low <- clean("base105low.RDATA","base105low")
base105high <- clean("base105high.RDATA","base105high")
base105scenarios <- bindScenarios(base105low, base105high)
# Other experimental scenarios
savePELPSToRdata("pelps2016noexchg", compressed = "zip")
savePELPSToRdata("pelps2016pmgnonews", compressed = "zip")
savePELPSToRdata("pelps2016dolsnonewsnoswd", compressed = "zip")
savePELPSToRdata("pelps2016simangunsong", compressed = "zip")
base2016noexchg <- clean("pelps2016noexchg.RDATA", "base2016noexchg")
pmgnonews <- clean("pelps2016pmgnonews.RDATA", "pmgnonews")
dolsnonewsnoswd <- clean("pelps2016dolsnonewsnoswd.RDATA",
"dolsnonewsnoswd")
simangunsong <- clean("pelps2016simangunsong.RDATA", "simangunsong")
message("Add or remove scenarios from this list for
future storage and visualisation as needed.")
scenarios <- Reduce(bindScenarios,
list(base2016noexchg, pmgnonews,
dolsnonewsnoswd, simangunsong))
# Write scenarios to a rds file for later use
saveRDS(scenarios, "enddata/base_dols_pmg2016.rds")
scenarios <- readRDS("enddata/base_dols_pmg2016.rds")
# scenarios <- readRDS("enddata/baseandpc2016.rds")
if(FALSE){
# Load TTIP high and low scenarios to compare
load("enddata/GFPM_Output_TTIP_with_sensitivity.RDATA")
scenarios <- allScenarios
# Load the re-run version of these scenarios
scenarios <- base105scenarios
}
GFPM simulation results
Plot consumption
# Plot demand for all products
p <- ggplot(data = NULL) +
aes(x = Period, y = Volume,
colour = GFPM_REG, linetype = Scenario) +
ggtitle("Demand") +
facet_wrap(~Product, scales="free_y") +
theme(legend.position = "bottom")
p + geom_line(data = filter(scenarios$aggregates,
Element == "Demand"))
# Plot demand for selected products only
selectedproducts4plot <- selectedproducts
# selectedproducts4plot <- "OthPaper"
p + geom_line(data = filter(scenarios$aggregates,
Element == "Demand" &
Product %in% selectedproducts4plot))
# Europe only
scenarios$aggregates %>%
filter(Element == "Demand" &
GFPM_REG == "Europe" &
Product %in%
selectedproducts4plot) %>%
mutate(Year = Period * 5 + 2012) %>%
ggplot() +
aes(x = Year, y = Volume/1e3,
color = Scenario, linetype = Scenario) +
# ggtitle("GFPM Demand scenarios for Europe") +
ylab("Volume in million T or m3") +
facet_wrap(~Product, scales="free_y") +
theme(legend.position = "bottom") +
geom_line(data = ) +
ylim(0,NA)
# # Europe only without Simangunsong
# p + geom_line(data = filter(scenarios$aggregates,
# Element == "Demand" &
# GFPM_REG == "Europe" &
# Product %in%
# selectedproducts4plot &
# Scenario != "simangunsong"),
# mapping = aes(color = Scenario)) +
# ylim(0,NA)
#
EUcountries <- filter(countrycodes,EU27)$Country
eugfpm <- scenarios$entity %>%
filter(Country %in% EUcountries & Element == "Demand" &
Product %in% selectedproducts) %>%
# Rename Scenarios so that they match the elasticities table in the article
rename(Scenario2replace = Scenario) %>%
left_join(data_frame(Scenario2replace = c("base2016noexchg", "dolsnonewsnoswd",
"pmgnonews", "simangunsong"),
Scenario = c("Buongiorno 2015", "DOLS", "PMG", "Simangunsong 2001")),
by="Scenario2replace") %>%
# Rename Items so that they matche the elasticities table
left_join(productnames, by = "Product") %>%
group_by(Scenario, item, Period) %>%
summarise(Volume = sum(Volume)) %>%
mutate(Year = (Period-1) * 5 + 2012)
eugfpm %>%
ggplot() +
aes(x = Year, y = Volume/1e3,
color = Scenario, linetype = Scenario) +
# ggtitle("GFPM Demand scenarios for Europe") +
ylab(expression(paste("Volume in million ", m^{3},"| Volume in million T"))) +
facet_wrap(~item, scales="free_y") +
theme(legend.position = "bottom") +
geom_line(data = ) +
ylim(0,NA)
p <- ggplot(data = NULL) +
aes(x = Period, y = Volume,
colour = Country, linetype = Scenario) +
ggtitle("Demand") +
facet_wrap(~Product, scales="free_y") +
theme(legend.position = "bottom")
p + geom_line(data=filter(scenarios$entity,
Element == "Demand" &
GFPM_REG == "Europe" &
Product %in%
selectedproducts4plot))
Differences between base, DOLS and PMG scenarios
# Are values exactly the same?
scendiff <- scenarios$entity %>%
spread(Scenario,Volume) %>%
mutate(diffdols = base2016noexchg - dolsnonewsnoswd,
diffpmg = base2016noexchg - pmgnonews) %>%
arrange(diffdols)
# head(scendiff,100)
# tail(scendiff,100)
ggplot(data = filter(scendiff,
Country %in% c("France", "Germany",
"United Kingdom"),
Element == "Demand")) +
geom_line(mapping = aes(x = Period, y = diffdols, colour = Country, linetype = "dols")) +
geom_line(mapping = aes(x = Period, y = diffpmg, colour = Country, linetype = "pmg")) +
ggtitle("Demand") +
ylab("Difference with the base scenario") +
theme(legend.position = "bottom") +
facet_wrap(~Product, scales="free_y")
# Aggregated differences by region
scendiffagg <- scenarios$aggregates %>%
spread(Scenario, Volume) %>%
mutate(diffdols = base2016noexchg - dolsnonewsnoswd,
diffdolspercent = round(diffdols / base2016noexchg*100),
diffpmg = base2016noexchg - pmgnonews,
diffpmgpercent = round(diffpmg / base2016noexchg*100))
# scendiffagg %>% arrange(diff) %>% tail()
# Difference in the 5th period by product by region
scendiffagg %>%
filter(Period == 5 &
Product %in% selectedproducts &
Element == "Demand")
# # Same in percent
# scendiffagg %>%
# filter(Period == 5 &
# Product %in% selectedproducts &
# Element == "Demand")
scendiffagg %>%
filter(Element == "Demand" &
Period ==5 &
GFPM_REG == "Europe") %>%
select(Product, starts_with("diff")) %>%
kable()
# Plot
ggplot(data = filter(scendiffagg,
Element == "Demand")) +
aes(x = Period, colour = GFPM_REG) +
geom_line(aes(y = diffdols, linetype = "dols")) +
geom_line(aes(y = diffpmg, linetype = "pmg")) +
ggtitle("Demand") +
theme(legend.position = "bottom") + facet_wrap(~Product, scales="free_y")
Consumption table
Change in consumption in the EU in period 5
EUcountries <- filter(countrycodes,EU27)$Country
eu2032 <- scenarios$entity %>%
filter(Country %in% EUcountries &
Period == 5 &
Element == "Demand") %>%
group_by(Product, Scenario) %>%
summarise(Volume = round(sum(Volume)/1e3,3)) %>%
spread(Scenario, Volume) %>%
left_join(productnames, by="Product") %>%
filter(!is.na(item)) %>%
# Arrange items in the usual order
# mutate(item = factor(item, levels = productnames$item, ordered = TRUE)) %>%
ungroup() %>%
arrange(item) %>%
select(item, 2:5)
eu2032 %>%
kable()
eu2032 %>%
mutate(diffdols = (dolsnonewsnoswd - base2016noexchg)/base2016noexchg,
diffpmg = (pmgnonews -base2016noexchg)/base2016noexchg,
sdiffdols = (dolsnonewsnoswd - simangunsong)/simangunsong,
sdiffpmg = (pmgnonews -simangunsong)/simangunsong) %>%
select(item, 6:9)
paul4forest/GFPMoutput documentation built on May 24, 2019, 8:25 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(knitr) opts_knit$set(root.dir="../..") # file paths are relative to the root of the project directory library(GFPMoutput) library(dplyr) library(tidyr) library(ggplot2) # Attention do not load the tradeflow package here, # It's clean function masks this one. selectedproducts <- c("Newsprint", "PWPaper", "OthPaper", "Sawnwood", "ParticleB", "FiberB")
Blog
- Before July first try two different base files the default 2016 world file and base105 from previous work with Joseph.
- 20160706 try several ways to run the simulation It seems there is an error after the first period when the Excel files where modified in a French version of Excel. I changed the region and Language settings to English US. GFPM also gives also an error when there are positive price elasticities. I decided to only export those products which have negative price elasticities. The issue is that simulation results are exactly the same.
Load demand elasticities from the GFPM base scenario
worlddemand <- read.csv("rawdata/World2016demand.csv", stringsAsFactors = FALSE) %>% left_join(countrycodes[c("Country_Code","Country")]) %>% left_join(productcodes) worlddemand105 <- read.csv("rawdata/Worlddemand105.csv", stringsAsFactors = FALSE) %>% left_join(countrycodes[c("Country_Code","Country")]) %>% left_join(productcodes) # Country with zero elast worlddemand %>% filter(priceelast == 0) message("Dummy country zy can be ignored") # elastsummmary <- worlddemand %>% filter(priceelast != 0) %>% select(Product, priceelast, gdpelast) %>% distinct() elastsummmary %>% kable() elastsummmary %>% filter(Product %in% selectedproducts & priceelast !=-0.25) %>% dput()
Show products that have more than one elasticity
elastsummmary %>% filter(duplicated(Product)) %>% kable()
Show the group of countries which share a same elasticity of fuelwood demand.
fwdgroup <- worlddemand %>% filter(Product == "Fuelwood" & !is.na(Country)) %>% group_by(priceelast,gdpelast) %>% summarise(Country = paste(Country,collapse=",")) fwdgroup %>% kable()
plygroup <- worlddemand %>% filter(Product == "Plywood" & !is.na(Country)) %>% group_by(priceelast, gdpelast) %>% summarise(Country = paste(Country, collapse=",")) plygroup %>% kable()
newsgroup <- worlddemand %>% filter(Product == "Newsprint" & !is.na(Country)) %>% group_by(priceelast, gdpelast) %>% summarise(Country = paste(Country, collapse=",")) newsgroup %>% kable()
Update the elasticities table for the World.xls of GFPM
Use elasticities from the pannel cointegration (PC) article .
productnames <- data_frame(Product = c("Newsprint", "PWPaper", "OthPaper", "Sawnwood", "ParticleB", "FiberB", "Plywood"), item = c("Newsprint", "PrintingandWritingPaper", "OtherPaperandPaperboard", "SawnwoodConiferous", "ParticleBoard", "Fibreboard", "Plywood")) %>% mutate(item = factor(item, levels = item, ordered = TRUE)) # GDP, Price as in the article's table elast <- read.csv("/home/paul/R/forestproductsdemand/data-end/elastforGFPM.csv", stringsAsFactors = FALSE) %>% rename(gdpelastpc = lgdpreur, priceelastpc = lpricereur) %>% # Rename Items to GFPM Products left_join(productnames, by = "item") %>% select(-item) elast %>% kable() # Simangunsong product name missing before sawnwood # not important Total swd and total panel not useful anyway # Only plywood would be interesting. # Compare with elasticites currently in the GFPM elast %>% left_join(filter(elastsummmary, Product %in% elast$Product), by="Product") %>% select(Product, estimator, gdpelast, gdpelastpc, priceelast, priceelastpc) %>% kable()
# Chunk set to eval=FALSE because it's used only once before a simulation elast %>% filter(estimator == "dols") %>% adddemandelast(worlddemand105) %>% write.csv("rawdata/World105demanddols.csv") elast %>% filter(estimator == "pmg") %>% adddemandelast(worlddemand105) %>% write.csv("rawdata/World105demandpmg.csv") elast %>% filter(estimator == "dols") %>% adddemandelast(worlddemand) %>% write.csv("rawdata/World2016demanddols.csv") elast %>% filter(estimator == "pmg") %>% adddemandelast(worlddemand) %>% write.csv("rawdata/World2016demandpmg.csv") elast %>% filter(estimator == "simangunsong") %>% adddemandelast(worlddemand) %>% write.csv("rawdata/World2016demandsimangunsong.csv") if(FALSE){ # Only Sawnwood message("DOLS::Newsprint and Sawnwood have a positive price elasticity") message("PMG::Newsprint has a positive price elasticity") # DOLS elast without Newsprint without sawnwood elast %>% filter(estimator == "dols" & ! Product %in% c("Newsprint","Sawnwood")) %>% adddemandelast(worlddemand) %>% write.csv("rawdata/World2016demandnonewsnoswd.csv") # Only particle board elast %>% filter(estimator == "dols" & Product == "ParticleB") %>% adddemandelast(worlddemand) %>% write.csv("rawdata/World2016demanddolsParticleB.csv") # Only Othpaper elast %>% filter(estimator == "dols" & Product == "OthPaper") %>% adddemandelast(worlddemand) %>% write.csv("rawdata/World2016demanddolsOthPaper.csv") # Pmg elast without newsprint, # Avoid newsprint because it's positive price elasticity # seems to be creating an error elast %>% filter(estimator == "pmg" & Product != "Newsprint") %>% adddemandelast(worlddemand) %>% write.csv("rawdata/World2016demandpmgnonews.csv") } if (FALSE){ pcelast <- elast %>% filter(estimator == "dols") # Develop the adddemandelast function worlddemand2 <- worlddemand %>% left_join(pcelast, by="Product") # Replace those that have a new elasticity by the new value worlddemand2 <- worlddemand2 %>% filter(!is.na(gdpelastpc)) %>% mutate(gdpelast = gdpelastpc, priceelast = priceelastpc) %>% # Add back those that didn't have a new elasticity rbind(filter(worlddemand2,is.na(gdpelastpc))) %>% arrange(Country_Code, Product_Code) # Did the number of rows stay the same? stopifnot(identical(nrow(worlddemand),nrow(worlddemand2))) # Did the order of rows stay the same? stopifnot(identical(worlddemand[c("Country_Code","Product_Code","basedemand")],worlddemand2[c("Country_Code","Product_Code","basedemand")])) write.csv(worlddemand2, "rawdata/World2016demandPC.csv", row.names = FALSE) worlddemand2016dols <- elast %>% filter(estimator == "dols") %>% adddemandelast(worlddemand) }
GFPM issues
Positive price elasticities
Positive price elasticities cause the GFPM simulation to break after the first optimisation round.
- The simulation breaks when I use all PMG estimates including newsprint elasticitity of demand.
- The simulation runs find when I use PMG estimates excluding newsprint elasticity of demand
- Similar issue for DOLS estimatesof the newsprint and sawnwood elasticity of demand.
Therefore I used only negative price elasticities. Price elasticities should be negative according to the theory.
Remove Exogchange elasticities changes
Joseph Buongiorno:
The reason for the small difference that you got is that in the world file that you used, the elasticities changed over time. This is reflected in the ExogChange file where the elasticities differ in each period. It looks like you have not changed that. So, the results reflect only the effect of a change in elasticity in the first period, which is small. I suggest changing the elasticities in all periods in ExogChange so that they are the same as in the first period and equal to your revised elasticities.
I removed GDP and price elasticity from exogenous change. By setting the GDP elasticity (column R) to 0 with a double click. Double click copies a full column down untill the next empty cell. Then I copied this column of 0 to the price elasticity (column H). Otherwise a double click on the price elastcitity section would errase all content from column H, for all data types. WHich is not what I want. I wand to modify the D data types only.
Excel issues
Because I use a French/ English version of excel I had to change dots to commas in the file. I swidched windows regional and language settings to US so that I don't have to do this commas / dot replacement each time I load a csv file.
Then here is how to send updated elasticities to GFPM:
- open this csv in Excel, convert the text to columns if necessary using comma as a separator.
- paste the priceelast and gdpelast columns in the World.xls table
- replace points by commas because we are dealing with a french version of excel (Side note: commas are incompatible with csv, they would require every field to be quoted. Now I understand why read.csv2 accepts text separated by semi-colons)
Then run the simulation as an alternative scenario and go back to this document to load both the base and alternative scenarios below.
Load simulation data
I copied the file to an archive called "pelpsbase2016"
under ./rawdata.
This needs to be evaluated only once.
Once data has been processed and storred in enddata.
It can be read from there.
message("Chunk set to eval=FALSE as it needs to be evaluated only once to load data from GFPM after a simulation.") savePELPSToRdata("pelpsbase2016", compressed = "zip") savePELPSToRdata("pelpspc2016", compressed = "zip") savePELPSToRdata("pelpslow2016", compressed = "zip") base2016 <- clean("pelpsbase2016.RDATA","base2016") pc <- clean("pelpspc2016.RDATA","pc2016") low <- clean("pelpslow2016.RDATA","low2016") # Write to endata message("Use bindScenarios() to combine scenarios ", "in one list, then write the two scenarios in one rds file.") scenarios <- Reduce(bindScenarios, list(base, pc, low)) saveRDS(scenarios, "enddata/baseandpc2016.rds") # base105 scenarios savePELPSToRdata("base105low", compressed = "zip") savePELPSToRdata("base105high", compressed = "zip") base105low <- clean("base105low.RDATA","base105low") base105high <- clean("base105high.RDATA","base105high") base105scenarios <- bindScenarios(base105low, base105high) # Other experimental scenarios savePELPSToRdata("pelps2016noexchg", compressed = "zip") savePELPSToRdata("pelps2016pmgnonews", compressed = "zip") savePELPSToRdata("pelps2016dolsnonewsnoswd", compressed = "zip") savePELPSToRdata("pelps2016simangunsong", compressed = "zip") base2016noexchg <- clean("pelps2016noexchg.RDATA", "base2016noexchg") pmgnonews <- clean("pelps2016pmgnonews.RDATA", "pmgnonews") dolsnonewsnoswd <- clean("pelps2016dolsnonewsnoswd.RDATA", "dolsnonewsnoswd") simangunsong <- clean("pelps2016simangunsong.RDATA", "simangunsong") message("Add or remove scenarios from this list for future storage and visualisation as needed.") scenarios <- Reduce(bindScenarios, list(base2016noexchg, pmgnonews, dolsnonewsnoswd, simangunsong)) # Write scenarios to a rds file for later use saveRDS(scenarios, "enddata/base_dols_pmg2016.rds")
scenarios <- readRDS("enddata/base_dols_pmg2016.rds") # scenarios <- readRDS("enddata/baseandpc2016.rds") if(FALSE){ # Load TTIP high and low scenarios to compare load("enddata/GFPM_Output_TTIP_with_sensitivity.RDATA") scenarios <- allScenarios # Load the re-run version of these scenarios scenarios <- base105scenarios }
GFPM simulation results
Plot consumption
# Plot demand for all products p <- ggplot(data = NULL) + aes(x = Period, y = Volume, colour = GFPM_REG, linetype = Scenario) + ggtitle("Demand") + facet_wrap(~Product, scales="free_y") + theme(legend.position = "bottom") p + geom_line(data = filter(scenarios$aggregates, Element == "Demand")) # Plot demand for selected products only selectedproducts4plot <- selectedproducts # selectedproducts4plot <- "OthPaper" p + geom_line(data = filter(scenarios$aggregates, Element == "Demand" & Product %in% selectedproducts4plot)) # Europe only scenarios$aggregates %>% filter(Element == "Demand" & GFPM_REG == "Europe" & Product %in% selectedproducts4plot) %>% mutate(Year = Period * 5 + 2012) %>% ggplot() + aes(x = Year, y = Volume/1e3, color = Scenario, linetype = Scenario) + # ggtitle("GFPM Demand scenarios for Europe") + ylab("Volume in million T or m3") + facet_wrap(~Product, scales="free_y") + theme(legend.position = "bottom") + geom_line(data = ) + ylim(0,NA) # # Europe only without Simangunsong # p + geom_line(data = filter(scenarios$aggregates, # Element == "Demand" & # GFPM_REG == "Europe" & # Product %in% # selectedproducts4plot & # Scenario != "simangunsong"), # mapping = aes(color = Scenario)) + # ylim(0,NA) #
EUcountries <- filter(countrycodes,EU27)$Country eugfpm <- scenarios$entity %>% filter(Country %in% EUcountries & Element == "Demand" & Product %in% selectedproducts) %>% # Rename Scenarios so that they match the elasticities table in the article rename(Scenario2replace = Scenario) %>% left_join(data_frame(Scenario2replace = c("base2016noexchg", "dolsnonewsnoswd", "pmgnonews", "simangunsong"), Scenario = c("Buongiorno 2015", "DOLS", "PMG", "Simangunsong 2001")), by="Scenario2replace") %>% # Rename Items so that they matche the elasticities table left_join(productnames, by = "Product") %>% group_by(Scenario, item, Period) %>% summarise(Volume = sum(Volume)) %>% mutate(Year = (Period-1) * 5 + 2012) eugfpm %>% ggplot() + aes(x = Year, y = Volume/1e3, color = Scenario, linetype = Scenario) + # ggtitle("GFPM Demand scenarios for Europe") + ylab(expression(paste("Volume in million ", m^{3},"| Volume in million T"))) + facet_wrap(~item, scales="free_y") + theme(legend.position = "bottom") + geom_line(data = ) + ylim(0,NA)
p <- ggplot(data = NULL) + aes(x = Period, y = Volume, colour = Country, linetype = Scenario) + ggtitle("Demand") + facet_wrap(~Product, scales="free_y") + theme(legend.position = "bottom") p + geom_line(data=filter(scenarios$entity, Element == "Demand" & GFPM_REG == "Europe" & Product %in% selectedproducts4plot))
Differences between base, DOLS and PMG scenarios
# Are values exactly the same? scendiff <- scenarios$entity %>% spread(Scenario,Volume) %>% mutate(diffdols = base2016noexchg - dolsnonewsnoswd, diffpmg = base2016noexchg - pmgnonews) %>% arrange(diffdols) # head(scendiff,100) # tail(scendiff,100) ggplot(data = filter(scendiff, Country %in% c("France", "Germany", "United Kingdom"), Element == "Demand")) + geom_line(mapping = aes(x = Period, y = diffdols, colour = Country, linetype = "dols")) + geom_line(mapping = aes(x = Period, y = diffpmg, colour = Country, linetype = "pmg")) + ggtitle("Demand") + ylab("Difference with the base scenario") + theme(legend.position = "bottom") + facet_wrap(~Product, scales="free_y") # Aggregated differences by region scendiffagg <- scenarios$aggregates %>% spread(Scenario, Volume) %>% mutate(diffdols = base2016noexchg - dolsnonewsnoswd, diffdolspercent = round(diffdols / base2016noexchg*100), diffpmg = base2016noexchg - pmgnonews, diffpmgpercent = round(diffpmg / base2016noexchg*100)) # scendiffagg %>% arrange(diff) %>% tail() # Difference in the 5th period by product by region scendiffagg %>% filter(Period == 5 & Product %in% selectedproducts & Element == "Demand") # # Same in percent # scendiffagg %>% # filter(Period == 5 & # Product %in% selectedproducts & # Element == "Demand") scendiffagg %>% filter(Element == "Demand" & Period ==5 & GFPM_REG == "Europe") %>% select(Product, starts_with("diff")) %>% kable() # Plot ggplot(data = filter(scendiffagg, Element == "Demand")) + aes(x = Period, colour = GFPM_REG) + geom_line(aes(y = diffdols, linetype = "dols")) + geom_line(aes(y = diffpmg, linetype = "pmg")) + ggtitle("Demand") + theme(legend.position = "bottom") + facet_wrap(~Product, scales="free_y")
Consumption table
Change in consumption in the EU in period 5
EUcountries <- filter(countrycodes,EU27)$Country eu2032 <- scenarios$entity %>% filter(Country %in% EUcountries & Period == 5 & Element == "Demand") %>% group_by(Product, Scenario) %>% summarise(Volume = round(sum(Volume)/1e3,3)) %>% spread(Scenario, Volume) %>% left_join(productnames, by="Product") %>% filter(!is.na(item)) %>% # Arrange items in the usual order # mutate(item = factor(item, levels = productnames$item, ordered = TRUE)) %>% ungroup() %>% arrange(item) %>% select(item, 2:5) eu2032 %>% kable() eu2032 %>% mutate(diffdols = (dolsnonewsnoswd - base2016noexchg)/base2016noexchg, diffpmg = (pmgnonews -base2016noexchg)/base2016noexchg, sdiffdols = (dolsnonewsnoswd - simangunsong)/simangunsong, sdiffpmg = (pmgnonews -simangunsong)/simangunsong) %>% select(item, 6:9)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.