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")
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()
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()
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) }
Positive price elasticities cause the GFPM simulation to break after the first optimisation round.
Therefore I used only negative price elasticities. Price elasticities should be negative according to the theory.
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.
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:
Then run the simulation as an alternative scenario and go back to this document to load both the base and alternative scenarios below.
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 }
# 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))
# 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")
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.