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

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.

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:

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.