####################################################################
## author:    Michael L. Wicki
## contact:   wimi@ethz.ch, ETH Zurich
## file name: mobility_conjoint.R
## Context:   Data Analysis, Robustness Checks
## started:   2018-12-17
####################################################################

#Copyright (C) 2020 Michael Wicki

#This R-Script is free software: you can redistribute it and/or modify
#it under the terms of the GNU General Public License as published by
#the Free Software Foundation, either version 3 of the License, or
#any later version.

#This program is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#GNU General Public License for more details.

#You should have received a copy of the GNU General Public License
#along with this program.  If not, see <http://www.gnu.org/licenses/>.

rm(list = ls())
getwd()
setwd("\\\\gess-fs.d.ethz.ch\\home$\\wimi\\Documents\\Sufficiency\\Conjoint")
set.seed(42)
.libPaths( "\\\\gess-fs.d.ethz.ch\\home$\\wimi\\Documents\\R\\win-library\\3.3" )
#.libPaths( "\\\\gess-fs.d.ethz.ch\\home$\\wimi\\R3UserLibs" )



################################ load packages ################################ 
#install.packages("Rmisc")


library(Matrix)
library(ggthemes)
library(cjoint)
library(ggplot2)
library(effects)
library(lme4)
library(data.table)
library(dplyr)
library(nFactors)
library(texreg)
library(tidyr)
library(forcats)
library(Rmisc)
library(cowplot)
library(scales)


################################ load dataset ################################ 

load("conjoint_design.RData")


################## Choice by country, highest income classes ###############
################## DE ###############
df_conj_DE <- subset(df_conj, country == "DE" & incgrp=="4th Quintile" | country == "DE" & incgrp=="5th Quintile")
c_base_DE <- amce(choice ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab +  
                    attrib5_lab +  attrib4_lab +  attrib6_lab, data = df_conj_DE, 
                  cluster=TRUE, respondent.id="id", 
                  design=design, baselines = baselines_base)


plot(c_base_DE)


c_attribute1 <- c(0, c_base_DE$estimates$attrib1lab["AMCE","attrib1lab2Increasingpricesby15"], c_base_DE$estimates$attrib1lab["AMCE","attrib1lab1Increasingpricesby30"])
s_attribute1 <- c(0, c_base_DE$estimates$attrib1lab["Std. Error","attrib1lab2Increasingpricesby15"], c_base_DE$estimates$attrib1lab["Std. Error","attrib1lab1Increasingpricesby30"])
c_attribute2 <- c(0, c_base_DE$estimates$attrib2lab["AMCE","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], c_base_DE$estimates$attrib2lab["AMCE","attrib2lab3Publicprogramsforlowincomehouseholds"],c_base_DE$estimates$attrib2lab["AMCE","attrib2lab4Reduceincometaxes"], c_base_DE$estimates$attrib2lab["AMCE","attrib2lab5Notaxrevenues"])
s_attribute2 <- c(0, c_base_DE$estimates$attrib2lab["Std. Error","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], c_base_DE$estimates$attrib2lab["Std. Error","attrib2lab3Publicprogramsforlowincomehouseholds"],c_base_DE$estimates$attrib2lab["Std. Error","attrib2lab4Reduceincometaxes"],c_base_DE$estimates$attrib2lab["Std. Error","attrib2lab5Notaxrevenues"])
c_attribute7 <- c(0, c_base_DE$estimates$attrib7lab["AMCE","attrib7lab2Halvingsubsidies"], c_base_DE$estimates$attrib7lab["AMCE","attrib7lab1Eliminatingsubsidies"])
s_attribute7 <- c(0, c_base_DE$estimates$attrib7lab["Std. Error","attrib7lab2Halvingsubsidies"], c_base_DE$estimates$attrib7lab["Std. Error","attrib7lab1Eliminatingsubsidies"])
c_attribute3 <- c(0, c_base_DE$estimates$attrib3lab["AMCE","attrib3lab41dayperweek"], c_base_DE$estimates$attrib3lab["AMCE","attrib3lab33daysperweek"], c_base_DE$estimates$attrib3lab["AMCE","attrib3lab25daysperweek"], c_base_DE$estimates$attrib3lab["AMCE","attrib3lab1Alwaysrestricted"])
s_attribute3 <- c(0, c_base_DE$estimates$attrib3lab["Std. Error","attrib3lab41dayperweek"], c_base_DE$estimates$attrib3lab["Std. Error","attrib3lab33daysperweek"], c_base_DE$estimates$attrib3lab["Std. Error","attrib3lab25daysperweek"], c_base_DE$estimates$attrib3lab["Std. Error","attrib3lab1Alwaysrestricted"])
c_attribute5 <- c(0, c_base_DE$estimates$attrib5lab["AMCE","attrib5lab2Atleast40milespergallon"], c_base_DE$estimates$attrib5lab["AMCE","attrib5lab1Atleast55milespergallon"])
s_attribute5 <- c(0, c_base_DE$estimates$attrib5lab["Std. Error","attrib5lab2Atleast40milespergallon"], c_base_DE$estimates$attrib5lab["Std. Error","attrib5lab1Atleast55milespergallon"])
c_attribute4 <- c(0, c_base_DE$estimates$attrib4lab["AMCE","attrib4lab2Reducingpricesby15"], c_base_DE$estimates$attrib4lab["AMCE","attrib4lab1Reducingpricesby30"])
s_attribute4 <- c(0, c_base_DE$estimates$attrib4lab["Std. Error","attrib4lab2Reducingpricesby15"], c_base_DE$estimates$attrib4lab["Std. Error","attrib4lab1Reducingpricesby30"])
c_attribute6 <- c(0, c_base_DE$estimates$attrib6lab["AMCE","attrib6lab2Occasionalcampaigns"], c_base_DE$estimates$attrib6lab["AMCE","attrib6lab1Frequentcampaigns"])
s_attribute6 <- c(0, c_base_DE$estimates$attrib6lab["Std. Error","attrib6lab2Occasionalcampaigns"], c_base_DE$estimates$attrib6lab["Std. Error","attrib6lab1Frequentcampaigns"])


c <- c(c_attribute1, c_attribute2, c_attribute7, c_attribute3, c_attribute5, c_attribute4, c_attribute6)
s <- c(s_attribute1, s_attribute2, s_attribute7, s_attribute3, s_attribute5, s_attribute4, s_attribute6)

df_c_base_DE <- data.frame(c,s)
df_c_base_DE$att <- c(c_base_DE$attributes$attrib1lab[1], c_base_DE$attributes$attrib1lab[2:3], 
                        c_base_DE$attributes$attrib2lab[1], c_base_DE$attributes$attrib2lab[2:5],
                        c_base_DE$attributes$attrib7lab[3], c_base_DE$attributes$attrib7lab[2], c_base_DE$attributes$attrib7lab[1],
                        c_base_DE$attributes$attrib3lab[5], c_base_DE$attributes$attrib3lab[4], c_base_DE$attributes$attrib3lab[3], c_base_DE$attributes$attrib3lab[2], c_base_DE$attributes$attrib3lab[1], 
                        c_base_DE$attributes$attrib5lab[3], c_base_DE$attributes$attrib5lab[2], c_base_DE$attributes$attrib5lab[1],
                        c_base_DE$attributes$attrib4lab[3], c_base_DE$attributes$attrib4lab[2], c_base_DE$attributes$attrib4lab[1], 
                        c_base_DE$attributes$attrib6lab[3], c_base_DE$attributes$attrib6lab[2], c_base_DE$attributes$attrib6lab[1])
df_c_base_DE$att <- factor(df_c_base_DE$att, levels = c(rev(c("3Nonewtax", "2Increasingpricesby15", "1Increasingpricesby30",
                                                                  "2Generalgovernmentbudget", "1Publicenvironmentalandclimateprotectionprograms", "3Publicprogramsforlowincomehouseholds", "4Reduceincometaxes", "5Notaxrevenues",
                                                                  "3Keepingsubsidiesatcurrentlevel", "2Halvingsubsidies", "1Eliminatingsubsidies",
                                                                  "5Neverrestricted", "41dayperweek", "33daysperweek", "25daysperweek", "1Alwaysrestricted", 
                                                                  "3Nonewrequirements", "2Atleast40milespergallon", "1Atleast55milespergallon",
                                                                  "3Noincreaseofsupport", "2Reducingpricesby15", "1Reducingpricesby30",
                                                                  "3Nocampaigns", "2Occasionalcampaigns", "1Frequentcampaigns"))),
                             labels = c(rev(c("Baseline: No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                              "Baseline: General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                              "Baseline: Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                              "Baseline: Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                              "Baseline: No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                              "Baseline: No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                              "Baseline: No campaigns", "Occasional campaigns", "Frequent campaigns"))))

df_c_base_DE$policy <- ifelse(df_c_base_DE$att == "Baseline: No new tax" |  df_c_base_DE$att == "Increasing prices by 15%" |  df_c_base_DE$att == "Increasing prices by 30%", "1", ifelse(
  df_c_base_DE$att == "Baseline: General government budget" | df_c_base_DE$att == "Public environmental and climate protection programs" | df_c_base_DE$att == "Public programs for low-income households" | df_c_base_DE$att == "Reduce income taxes" | df_c_base_DE$att == "No tax revenues", "2", ifelse(
    df_c_base_DE$att == "Baseline: Keeping subsidies at current level" | df_c_base_DE$att == "Halving subsidies" | df_c_base_DE$att == "Eliminating subsidies", "3", ifelse(
      df_c_base_DE$att == "Baseline: Never restricted" | df_c_base_DE$att == "1 day per week" | df_c_base_DE$att == "3 days per week" | df_c_base_DE$att == "5 days per week" | df_c_base_DE$att == "Always restricted", "4", ifelse(
        df_c_base_DE$att == "Baseline: No new requirements" | df_c_base_DE$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_c_base_DE$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          df_c_base_DE$att == "Baseline: No campaigns" | df_c_base_DE$att == "Occasional campaigns" | df_c_base_DE$att == "Frequent campaigns", "6", ifelse(
            df_c_base_DE$att == "Baseline: No increase of support" | df_c_base_DE$att == "Reducing prices by 15%" | df_c_base_DE$att == "Reducing prices by 30%", "7",
            "8")))))))

df_c_base_DE$policy = factor(df_c_base_DE$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))

p_c_base_DE <- ggplot(df_c_base_DE, aes(x=att, y = c)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s)) +
  theme_bw() + 
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") 
p_c_base_DE

################## CN ###############

df_conj_CN <- subset(df_conj, country == "CN" & incgrp=="4th Quintile" | country == "CN" & incgrp=="5th Quintile")
c_base_CN <- amce(choice ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab +  
                    attrib5_lab +  attrib4_lab +  attrib6_lab, data = df_conj_CN, 
                  cluster=TRUE, respondent.id="id", 
                  design=design, baselines = baselines_base)


plot(c_base_CN)


c_attribute1 <- c(0, c_base_CN$estimates$attrib1lab["AMCE","attrib1lab2Increasingpricesby15"], c_base_CN$estimates$attrib1lab["AMCE","attrib1lab1Increasingpricesby30"])
s_attribute1 <- c(0, c_base_CN$estimates$attrib1lab["Std. Error","attrib1lab2Increasingpricesby15"], c_base_CN$estimates$attrib1lab["Std. Error","attrib1lab1Increasingpricesby30"])
c_attribute2 <- c(0, c_base_CN$estimates$attrib2lab["AMCE","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], c_base_CN$estimates$attrib2lab["AMCE","attrib2lab3Publicprogramsforlowincomehouseholds"],c_base_CN$estimates$attrib2lab["AMCE","attrib2lab4Reduceincometaxes"], c_base_CN$estimates$attrib2lab["AMCE","attrib2lab5Notaxrevenues"])
s_attribute2 <- c(0, c_base_CN$estimates$attrib2lab["Std. Error","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], c_base_CN$estimates$attrib2lab["Std. Error","attrib2lab3Publicprogramsforlowincomehouseholds"],c_base_CN$estimates$attrib2lab["Std. Error","attrib2lab4Reduceincometaxes"],c_base_CN$estimates$attrib2lab["Std. Error","attrib2lab5Notaxrevenues"])
c_attribute7 <- c(0, c_base_CN$estimates$attrib7lab["AMCE","attrib7lab2Halvingsubsidies"], c_base_CN$estimates$attrib7lab["AMCE","attrib7lab1Eliminatingsubsidies"])
s_attribute7 <- c(0, c_base_CN$estimates$attrib7lab["Std. Error","attrib7lab2Halvingsubsidies"], c_base_CN$estimates$attrib7lab["Std. Error","attrib7lab1Eliminatingsubsidies"])
c_attribute3 <- c(0, c_base_CN$estimates$attrib3lab["AMCE","attrib3lab41dayperweek"], c_base_CN$estimates$attrib3lab["AMCE","attrib3lab33daysperweek"], c_base_CN$estimates$attrib3lab["AMCE","attrib3lab25daysperweek"], c_base_CN$estimates$attrib3lab["AMCE","attrib3lab1Alwaysrestricted"])
s_attribute3 <- c(0, c_base_CN$estimates$attrib3lab["Std. Error","attrib3lab41dayperweek"], c_base_CN$estimates$attrib3lab["Std. Error","attrib3lab33daysperweek"], c_base_CN$estimates$attrib3lab["Std. Error","attrib3lab25daysperweek"], c_base_CN$estimates$attrib3lab["Std. Error","attrib3lab1Alwaysrestricted"])
c_attribute5 <- c(0, c_base_CN$estimates$attrib5lab["AMCE","attrib5lab2Atleast40milespergallon"], c_base_CN$estimates$attrib5lab["AMCE","attrib5lab1Atleast55milespergallon"])
s_attribute5 <- c(0, c_base_CN$estimates$attrib5lab["Std. Error","attrib5lab2Atleast40milespergallon"], c_base_CN$estimates$attrib5lab["Std. Error","attrib5lab1Atleast55milespergallon"])
c_attribute4 <- c(0, c_base_CN$estimates$attrib4lab["AMCE","attrib4lab2Reducingpricesby15"], c_base_CN$estimates$attrib4lab["AMCE","attrib4lab1Reducingpricesby30"])
s_attribute4 <- c(0, c_base_CN$estimates$attrib4lab["Std. Error","attrib4lab2Reducingpricesby15"], c_base_CN$estimates$attrib4lab["Std. Error","attrib4lab1Reducingpricesby30"])
c_attribute6 <- c(0, c_base_CN$estimates$attrib6lab["AMCE","attrib6lab2Occasionalcampaigns"], c_base_CN$estimates$attrib6lab["AMCE","attrib6lab1Frequentcampaigns"])
s_attribute6 <- c(0, c_base_CN$estimates$attrib6lab["Std. Error","attrib6lab2Occasionalcampaigns"], c_base_CN$estimates$attrib6lab["Std. Error","attrib6lab1Frequentcampaigns"])


c <- c(c_attribute1, c_attribute2, c_attribute7, c_attribute3, c_attribute5, c_attribute4, c_attribute6)
s <- c(s_attribute1, s_attribute2, s_attribute7, s_attribute3, s_attribute5, s_attribute4, s_attribute6)

df_c_base_CN <- data.frame(c,s)
df_c_base_CN$att <- c(c_base_CN$attributes$attrib1lab[1], c_base_CN$attributes$attrib1lab[2:3], 
                      c_base_CN$attributes$attrib2lab[1], c_base_CN$attributes$attrib2lab[2:5],
                      c_base_CN$attributes$attrib7lab[3], c_base_CN$attributes$attrib7lab[2], c_base_CN$attributes$attrib7lab[1],
                      c_base_CN$attributes$attrib3lab[5], c_base_CN$attributes$attrib3lab[4], c_base_CN$attributes$attrib3lab[3], c_base_CN$attributes$attrib3lab[2], c_base_CN$attributes$attrib3lab[1], 
                      c_base_CN$attributes$attrib5lab[3], c_base_CN$attributes$attrib5lab[2], c_base_CN$attributes$attrib5lab[1],
                      c_base_CN$attributes$attrib4lab[3], c_base_CN$attributes$attrib4lab[2], c_base_CN$attributes$attrib4lab[1], 
                      c_base_CN$attributes$attrib6lab[3], c_base_CN$attributes$attrib6lab[2], c_base_CN$attributes$attrib6lab[1])
df_c_base_CN$att <- factor(df_c_base_CN$att, levels = c(rev(c("3Nonewtax", "2Increasingpricesby15", "1Increasingpricesby30",
                                                              "2Generalgovernmentbudget", "1Publicenvironmentalandclimateprotectionprograms", "3Publicprogramsforlowincomehouseholds", "4Reduceincometaxes", "5Notaxrevenues",
                                                              "3Keepingsubsidiesatcurrentlevel", "2Halvingsubsidies", "1Eliminatingsubsidies",
                                                              "5Neverrestricted", "41dayperweek", "33daysperweek", "25daysperweek", "1Alwaysrestricted", 
                                                              "3Nonewrequirements", "2Atleast40milespergallon", "1Atleast55milespergallon",
                                                              "3Noincreaseofsupport", "2Reducingpricesby15", "1Reducingpricesby30",
                                                              "3Nocampaigns", "2Occasionalcampaigns", "1Frequentcampaigns"))),
                           labels = c(rev(c("Baseline: No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                            "Baseline: General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                            "Baseline: Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                            "Baseline: Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                            "Baseline: No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                            "Baseline: No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                            "Baseline: No campaigns", "Occasional campaigns", "Frequent campaigns"))))

df_c_base_CN$policy <- ifelse(df_c_base_CN$att == "Baseline: No new tax" |  df_c_base_CN$att == "Increasing prices by 15%" |  df_c_base_CN$att == "Increasing prices by 30%", "1", ifelse(
  df_c_base_CN$att == "Baseline: General government budget" | df_c_base_CN$att == "Public environmental and climate protection programs" | df_c_base_CN$att == "Public programs for low-income households" | df_c_base_CN$att == "Reduce income taxes" | df_c_base_CN$att == "No tax revenues", "2", ifelse(
    df_c_base_CN$att == "Baseline: Keeping subsidies at current level" | df_c_base_CN$att == "Halving subsidies" | df_c_base_CN$att == "Eliminating subsidies", "3", ifelse(
      df_c_base_CN$att == "Baseline: Never restricted" | df_c_base_CN$att == "1 day per week" | df_c_base_CN$att == "3 days per week" | df_c_base_CN$att == "5 days per week" | df_c_base_CN$att == "Always restricted", "4", ifelse(
        df_c_base_CN$att == "Baseline: No new requirements" | df_c_base_CN$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_c_base_CN$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          df_c_base_CN$att == "Baseline: No campaigns" | df_c_base_CN$att == "Occasional campaigns" | df_c_base_CN$att == "Frequent campaigns", "6", ifelse(
            df_c_base_CN$att == "Baseline: No increase of support" | df_c_base_CN$att == "Reducing prices by 15%" | df_c_base_CN$att == "Reducing prices by 30%", "7",
            "8")))))))

df_c_base_CN$policy = factor(df_c_base_CN$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))


p_c_base_CN <- ggplot(df_c_base_CN, aes(x=att, y = c)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s)) +
  theme_bw() + 
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") 
p_c_base_CN

################## US ###############

df_conj_US <- subset(df_conj, country == "US" & incgrp=="4th Quintile" | country == "US" & incgrp=="5th Quintile")
c_base_US <- amce(choice ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab +  
                    attrib5_lab +  attrib4_lab +  attrib6_lab, data = df_conj_US, 
                  cluster=TRUE, respondent.id="id", 
                  design=design, baselines = baselines_base)


plot(c_base_US)


c_attribute1 <- c(0, c_base_US$estimates$attrib1lab["AMCE","attrib1lab2Increasingpricesby15"], c_base_US$estimates$attrib1lab["AMCE","attrib1lab1Increasingpricesby30"])
s_attribute1 <- c(0, c_base_US$estimates$attrib1lab["Std. Error","attrib1lab2Increasingpricesby15"], c_base_US$estimates$attrib1lab["Std. Error","attrib1lab1Increasingpricesby30"])
c_attribute2 <- c(0, c_base_US$estimates$attrib2lab["AMCE","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], c_base_US$estimates$attrib2lab["AMCE","attrib2lab3Publicprogramsforlowincomehouseholds"],c_base_US$estimates$attrib2lab["AMCE","attrib2lab4Reduceincometaxes"], c_base_US$estimates$attrib2lab["AMCE","attrib2lab5Notaxrevenues"])
s_attribute2 <- c(0, c_base_US$estimates$attrib2lab["Std. Error","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], c_base_US$estimates$attrib2lab["Std. Error","attrib2lab3Publicprogramsforlowincomehouseholds"],c_base_US$estimates$attrib2lab["Std. Error","attrib2lab4Reduceincometaxes"],c_base_US$estimates$attrib2lab["Std. Error","attrib2lab5Notaxrevenues"])
c_attribute7 <- c(0, c_base_US$estimates$attrib7lab["AMCE","attrib7lab2Halvingsubsidies"], c_base_US$estimates$attrib7lab["AMCE","attrib7lab1Eliminatingsubsidies"])
s_attribute7 <- c(0, c_base_US$estimates$attrib7lab["Std. Error","attrib7lab2Halvingsubsidies"], c_base_US$estimates$attrib7lab["Std. Error","attrib7lab1Eliminatingsubsidies"])
c_attribute3 <- c(0, c_base_US$estimates$attrib3lab["AMCE","attrib3lab41dayperweek"], c_base_US$estimates$attrib3lab["AMCE","attrib3lab33daysperweek"], c_base_US$estimates$attrib3lab["AMCE","attrib3lab25daysperweek"], c_base_US$estimates$attrib3lab["AMCE","attrib3lab1Alwaysrestricted"])
s_attribute3 <- c(0, c_base_US$estimates$attrib3lab["Std. Error","attrib3lab41dayperweek"], c_base_US$estimates$attrib3lab["Std. Error","attrib3lab33daysperweek"], c_base_US$estimates$attrib3lab["Std. Error","attrib3lab25daysperweek"], c_base_US$estimates$attrib3lab["Std. Error","attrib3lab1Alwaysrestricted"])
c_attribute5 <- c(0, c_base_US$estimates$attrib5lab["AMCE","attrib5lab2Atleast40milespergallon"], c_base_US$estimates$attrib5lab["AMCE","attrib5lab1Atleast55milespergallon"])
s_attribute5 <- c(0, c_base_US$estimates$attrib5lab["Std. Error","attrib5lab2Atleast40milespergallon"], c_base_US$estimates$attrib5lab["Std. Error","attrib5lab1Atleast55milespergallon"])
c_attribute4 <- c(0, c_base_US$estimates$attrib4lab["AMCE","attrib4lab2Reducingpricesby15"], c_base_US$estimates$attrib4lab["AMCE","attrib4lab1Reducingpricesby30"])
s_attribute4 <- c(0, c_base_US$estimates$attrib4lab["Std. Error","attrib4lab2Reducingpricesby15"], c_base_US$estimates$attrib4lab["Std. Error","attrib4lab1Reducingpricesby30"])
c_attribute6 <- c(0, c_base_US$estimates$attrib6lab["AMCE","attrib6lab2Occasionalcampaigns"], c_base_US$estimates$attrib6lab["AMCE","attrib6lab1Frequentcampaigns"])
s_attribute6 <- c(0, c_base_US$estimates$attrib6lab["Std. Error","attrib6lab2Occasionalcampaigns"], c_base_US$estimates$attrib6lab["Std. Error","attrib6lab1Frequentcampaigns"])


c <- c(c_attribute1, c_attribute2, c_attribute7, c_attribute3, c_attribute5, c_attribute4, c_attribute6)
s <- c(s_attribute1, s_attribute2, s_attribute7, s_attribute3, s_attribute5, s_attribute4, s_attribute6)

df_c_base_US <- data.frame(c,s)
df_c_base_US$att <- c(c_base_US$attributes$attrib1lab[1], c_base_US$attributes$attrib1lab[2:3], 
                      c_base_US$attributes$attrib2lab[1], c_base_US$attributes$attrib2lab[2:5],
                      c_base_US$attributes$attrib7lab[3], c_base_US$attributes$attrib7lab[2], c_base_US$attributes$attrib7lab[1],
                      c_base_US$attributes$attrib3lab[5], c_base_US$attributes$attrib3lab[4], c_base_US$attributes$attrib3lab[3], c_base_US$attributes$attrib3lab[2], c_base_US$attributes$attrib3lab[1], 
                      c_base_US$attributes$attrib5lab[3], c_base_US$attributes$attrib5lab[2], c_base_US$attributes$attrib5lab[1],
                      c_base_US$attributes$attrib4lab[3], c_base_US$attributes$attrib4lab[2], c_base_US$attributes$attrib4lab[1], 
                      c_base_US$attributes$attrib6lab[3], c_base_US$attributes$attrib6lab[2], c_base_US$attributes$attrib6lab[1])
df_c_base_US$att <- factor(df_c_base_US$att, levels = c(rev(c("3Nonewtax", "2Increasingpricesby15", "1Increasingpricesby30",
                                                              "2Generalgovernmentbudget", "1Publicenvironmentalandclimateprotectionprograms", "3Publicprogramsforlowincomehouseholds", "4Reduceincometaxes", "5Notaxrevenues",
                                                              "3Keepingsubsidiesatcurrentlevel", "2Halvingsubsidies", "1Eliminatingsubsidies",
                                                              "5Neverrestricted", "41dayperweek", "33daysperweek", "25daysperweek", "1Alwaysrestricted", 
                                                              "3Nonewrequirements", "2Atleast40milespergallon", "1Atleast55milespergallon",
                                                              "3Noincreaseofsupport", "2Reducingpricesby15", "1Reducingpricesby30",
                                                              "3Nocampaigns", "2Occasionalcampaigns", "1Frequentcampaigns"))),
                           labels = c(rev(c("Baseline: No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                            "Baseline: General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                            "Baseline: Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                            "Baseline: Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                            "Baseline: No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                            "Baseline: No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                            "Baseline: No campaigns", "Occasional campaigns", "Frequent campaigns"))))

df_c_base_US$policy <- ifelse(df_c_base_US$att == "Baseline: No new tax" |  df_c_base_US$att == "Increasing prices by 15%" |  df_c_base_US$att == "Increasing prices by 30%", "1", ifelse(
  df_c_base_US$att == "Baseline: General government budget" | df_c_base_US$att == "Public environmental and climate protection programs" | df_c_base_US$att == "Public programs for low-income households" | df_c_base_US$att == "Reduce income taxes" | df_c_base_US$att == "No tax revenues", "2", ifelse(
    df_c_base_US$att == "Baseline: Keeping subsidies at current level" | df_c_base_US$att == "Halving subsidies" | df_c_base_US$att == "Eliminating subsidies", "3", ifelse(
      df_c_base_US$att == "Baseline: Never restricted" | df_c_base_US$att == "1 day per week" | df_c_base_US$att == "3 days per week" | df_c_base_US$att == "5 days per week" | df_c_base_US$att == "Always restricted", "4", ifelse(
        df_c_base_US$att == "Baseline: No new requirements" | df_c_base_US$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_c_base_US$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          df_c_base_US$att == "Baseline: No campaigns" | df_c_base_US$att == "Occasional campaigns" | df_c_base_US$att == "Frequent campaigns", "6", ifelse(
            df_c_base_US$att == "Baseline: No increase of support" | df_c_base_US$att == "Reducing prices by 15%" | df_c_base_US$att == "Reducing prices by 30%", "7",
            "8")))))))

df_c_base_US$policy = factor(df_c_base_US$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))


p_c_base_US <- ggplot(df_c_base_US, aes(x=att, y = c)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s)) +
  theme_bw() + 
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") 
p_c_base_US

################## Rating by country, highest income classes only ###############
################## DE ###############
df_conj_DE <- subset(df_conj, country == "DE" & incgrp=="4th Quintile" | country == "DE" & incgrp=="5th Quintile")
r_base_DE <- amce(rate ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab +  
                    attrib5_lab +  attrib4_lab +  attrib6_lab, data = df_conj_DE, 
                  cluster=TRUE, respondent.id="id", 
                  design=design, baselines = baselines_base)


plot(r_base_DE)


c_attribute1 <- c(0, r_base_DE$estimates$attrib1lab["AMCE","attrib1lab2Increasingpricesby15"], r_base_DE$estimates$attrib1lab["AMCE","attrib1lab1Increasingpricesby30"])
s_attribute1 <- c(0, r_base_DE$estimates$attrib1lab["Std. Error","attrib1lab2Increasingpricesby15"], r_base_DE$estimates$attrib1lab["Std. Error","attrib1lab1Increasingpricesby30"])
c_attribute2 <- c(0, r_base_DE$estimates$attrib2lab["AMCE","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], r_base_DE$estimates$attrib2lab["AMCE","attrib2lab3Publicprogramsforlowincomehouseholds"],r_base_DE$estimates$attrib2lab["AMCE","attrib2lab4Reduceincometaxes"], r_base_DE$estimates$attrib2lab["AMCE","attrib2lab5Notaxrevenues"])
s_attribute2 <- c(0, r_base_DE$estimates$attrib2lab["Std. Error","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], r_base_DE$estimates$attrib2lab["Std. Error","attrib2lab3Publicprogramsforlowincomehouseholds"],r_base_DE$estimates$attrib2lab["Std. Error","attrib2lab4Reduceincometaxes"],r_base_DE$estimates$attrib2lab["Std. Error","attrib2lab5Notaxrevenues"])
c_attribute7 <- c(0, r_base_DE$estimates$attrib7lab["AMCE","attrib7lab2Halvingsubsidies"], r_base_DE$estimates$attrib7lab["AMCE","attrib7lab1Eliminatingsubsidies"])
s_attribute7 <- c(0, r_base_DE$estimates$attrib7lab["Std. Error","attrib7lab2Halvingsubsidies"], r_base_DE$estimates$attrib7lab["Std. Error","attrib7lab1Eliminatingsubsidies"])
c_attribute3 <- c(0, r_base_DE$estimates$attrib3lab["AMCE","attrib3lab41dayperweek"], r_base_DE$estimates$attrib3lab["AMCE","attrib3lab33daysperweek"], r_base_DE$estimates$attrib3lab["AMCE","attrib3lab25daysperweek"], r_base_DE$estimates$attrib3lab["AMCE","attrib3lab1Alwaysrestricted"])
s_attribute3 <- c(0, r_base_DE$estimates$attrib3lab["Std. Error","attrib3lab41dayperweek"], r_base_DE$estimates$attrib3lab["Std. Error","attrib3lab33daysperweek"], r_base_DE$estimates$attrib3lab["Std. Error","attrib3lab25daysperweek"], r_base_DE$estimates$attrib3lab["Std. Error","attrib3lab1Alwaysrestricted"])
c_attribute5 <- c(0, r_base_DE$estimates$attrib5lab["AMCE","attrib5lab2Atleast40milespergallon"], r_base_DE$estimates$attrib5lab["AMCE","attrib5lab1Atleast55milespergallon"])
s_attribute5 <- c(0, r_base_DE$estimates$attrib5lab["Std. Error","attrib5lab2Atleast40milespergallon"], r_base_DE$estimates$attrib5lab["Std. Error","attrib5lab1Atleast55milespergallon"])
c_attribute4 <- c(0, r_base_DE$estimates$attrib4lab["AMCE","attrib4lab2Reducingpricesby15"], r_base_DE$estimates$attrib4lab["AMCE","attrib4lab1Reducingpricesby30"])
s_attribute4 <- c(0, r_base_DE$estimates$attrib4lab["Std. Error","attrib4lab2Reducingpricesby15"], r_base_DE$estimates$attrib4lab["Std. Error","attrib4lab1Reducingpricesby30"])
c_attribute6 <- c(0, r_base_DE$estimates$attrib6lab["AMCE","attrib6lab2Occasionalcampaigns"], r_base_DE$estimates$attrib6lab["AMCE","attrib6lab1Frequentcampaigns"])
s_attribute6 <- c(0, r_base_DE$estimates$attrib6lab["Std. Error","attrib6lab2Occasionalcampaigns"], r_base_DE$estimates$attrib6lab["Std. Error","attrib6lab1Frequentcampaigns"])


c <- c(c_attribute1, c_attribute2, c_attribute7, c_attribute3, c_attribute5, c_attribute4, c_attribute6)
s <- c(s_attribute1, s_attribute2, s_attribute7, s_attribute3, s_attribute5, s_attribute4, s_attribute6)

df_r_base_DE <- data.frame(c,s)
df_r_base_DE$att <- c(r_base_DE$attributes$attrib1lab[1], r_base_DE$attributes$attrib1lab[2:3], 
                      r_base_DE$attributes$attrib2lab[1], r_base_DE$attributes$attrib2lab[2:5],
                      r_base_DE$attributes$attrib7lab[3], r_base_DE$attributes$attrib7lab[2], r_base_DE$attributes$attrib7lab[1],
                      r_base_DE$attributes$attrib3lab[5], r_base_DE$attributes$attrib3lab[4], r_base_DE$attributes$attrib3lab[3], r_base_DE$attributes$attrib3lab[2], r_base_DE$attributes$attrib3lab[1], 
                      r_base_DE$attributes$attrib5lab[3], r_base_DE$attributes$attrib5lab[2], r_base_DE$attributes$attrib5lab[1],
                      r_base_DE$attributes$attrib4lab[3], r_base_DE$attributes$attrib4lab[2], r_base_DE$attributes$attrib4lab[1], 
                      r_base_DE$attributes$attrib6lab[3], r_base_DE$attributes$attrib6lab[2], r_base_DE$attributes$attrib6lab[1])
df_r_base_DE$att <- factor(df_r_base_DE$att, levels = c(rev(c("3Nonewtax", "2Increasingpricesby15", "1Increasingpricesby30",
                                                              "2Generalgovernmentbudget", "1Publicenvironmentalandclimateprotectionprograms", "3Publicprogramsforlowincomehouseholds", "4Reduceincometaxes", "5Notaxrevenues",
                                                              "3Keepingsubsidiesatcurrentlevel", "2Halvingsubsidies", "1Eliminatingsubsidies",
                                                              "5Neverrestricted", "41dayperweek", "33daysperweek", "25daysperweek", "1Alwaysrestricted", 
                                                              "3Nonewrequirements", "2Atleast40milespergallon", "1Atleast55milespergallon",
                                                              "3Noincreaseofsupport", "2Reducingpricesby15", "1Reducingpricesby30",
                                                              "3Nocampaigns", "2Occasionalcampaigns", "1Frequentcampaigns"))),
                           labels = c(rev(c("Baseline: No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                            "Baseline: General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                            "Baseline: Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                            "Baseline: Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                            "Baseline: No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                            "Baseline: No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                            "Baseline: No campaigns", "Occasional campaigns", "Frequent campaigns"))))

df_r_base_DE$policy <- ifelse(df_r_base_DE$att == "Baseline: No new tax" |  df_r_base_DE$att == "Increasing prices by 15%" |  df_r_base_DE$att == "Increasing prices by 30%", "1", ifelse(
  df_r_base_DE$att == "Baseline: General government budget" | df_r_base_DE$att == "Public environmental and climate protection programs" | df_r_base_DE$att == "Public programs for low-income households" | df_r_base_DE$att == "Reduce income taxes" | df_r_base_DE$att == "No tax revenues", "2", ifelse(
    df_r_base_DE$att == "Baseline: Keeping subsidies at current level" | df_r_base_DE$att == "Halving subsidies" | df_r_base_DE$att == "Eliminating subsidies", "3", ifelse(
      df_r_base_DE$att == "Baseline: Never restricted" | df_r_base_DE$att == "1 day per week" | df_r_base_DE$att == "3 days per week" | df_r_base_DE$att == "5 days per week" | df_r_base_DE$att == "Always restricted", "4", ifelse(
        df_r_base_DE$att == "Baseline: No new requirements" | df_r_base_DE$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_r_base_DE$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          df_r_base_DE$att == "Baseline: No campaigns" | df_r_base_DE$att == "Occasional campaigns" | df_r_base_DE$att == "Frequent campaigns", "6", ifelse(
            df_r_base_DE$att == "Baseline: No increase of support" | df_r_base_DE$att == "Reducing prices by 15%" | df_r_base_DE$att == "Reducing prices by 30%", "7",
            "8")))))))

df_r_base_DE$policy = factor(df_r_base_DE$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))

p_r_base_DE <- ggplot(df_r_base_DE, aes(x=att, y = c)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s)) +
  theme_bw() + 
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") 
p_r_base_DE

################## CN ###############

df_conj_CN <- subset(df_conj, country == "CN" & incgrp=="4th Quintile" | country == "CN" & incgrp=="5th Quintile")
r_base_CN <- amce(rate ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab +  
                    attrib5_lab +  attrib4_lab +  attrib6_lab, data = df_conj_CN, 
                  cluster=TRUE, respondent.id="id", 
                  design=design, baselines = baselines_base)


plot(r_base_CN)


c_attribute1 <- c(0, r_base_CN$estimates$attrib1lab["AMCE","attrib1lab2Increasingpricesby15"], r_base_CN$estimates$attrib1lab["AMCE","attrib1lab1Increasingpricesby30"])
s_attribute1 <- c(0, r_base_CN$estimates$attrib1lab["Std. Error","attrib1lab2Increasingpricesby15"], r_base_CN$estimates$attrib1lab["Std. Error","attrib1lab1Increasingpricesby30"])
c_attribute2 <- c(0, r_base_CN$estimates$attrib2lab["AMCE","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], r_base_CN$estimates$attrib2lab["AMCE","attrib2lab3Publicprogramsforlowincomehouseholds"],r_base_CN$estimates$attrib2lab["AMCE","attrib2lab4Reduceincometaxes"], r_base_CN$estimates$attrib2lab["AMCE","attrib2lab5Notaxrevenues"])
s_attribute2 <- c(0, r_base_CN$estimates$attrib2lab["Std. Error","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], r_base_CN$estimates$attrib2lab["Std. Error","attrib2lab3Publicprogramsforlowincomehouseholds"],r_base_CN$estimates$attrib2lab["Std. Error","attrib2lab4Reduceincometaxes"],r_base_CN$estimates$attrib2lab["Std. Error","attrib2lab5Notaxrevenues"])
c_attribute7 <- c(0, r_base_CN$estimates$attrib7lab["AMCE","attrib7lab2Halvingsubsidies"], r_base_CN$estimates$attrib7lab["AMCE","attrib7lab1Eliminatingsubsidies"])
s_attribute7 <- c(0, r_base_CN$estimates$attrib7lab["Std. Error","attrib7lab2Halvingsubsidies"], r_base_CN$estimates$attrib7lab["Std. Error","attrib7lab1Eliminatingsubsidies"])
c_attribute3 <- c(0, r_base_CN$estimates$attrib3lab["AMCE","attrib3lab41dayperweek"], r_base_CN$estimates$attrib3lab["AMCE","attrib3lab33daysperweek"], r_base_CN$estimates$attrib3lab["AMCE","attrib3lab25daysperweek"], r_base_CN$estimates$attrib3lab["AMCE","attrib3lab1Alwaysrestricted"])
s_attribute3 <- c(0, r_base_CN$estimates$attrib3lab["Std. Error","attrib3lab41dayperweek"], r_base_CN$estimates$attrib3lab["Std. Error","attrib3lab33daysperweek"], r_base_CN$estimates$attrib3lab["Std. Error","attrib3lab25daysperweek"], r_base_CN$estimates$attrib3lab["Std. Error","attrib3lab1Alwaysrestricted"])
c_attribute5 <- c(0, r_base_CN$estimates$attrib5lab["AMCE","attrib5lab2Atleast40milespergallon"], r_base_CN$estimates$attrib5lab["AMCE","attrib5lab1Atleast55milespergallon"])
s_attribute5 <- c(0, r_base_CN$estimates$attrib5lab["Std. Error","attrib5lab2Atleast40milespergallon"], r_base_CN$estimates$attrib5lab["Std. Error","attrib5lab1Atleast55milespergallon"])
c_attribute4 <- c(0, r_base_CN$estimates$attrib4lab["AMCE","attrib4lab2Reducingpricesby15"], r_base_CN$estimates$attrib4lab["AMCE","attrib4lab1Reducingpricesby30"])
s_attribute4 <- c(0, r_base_CN$estimates$attrib4lab["Std. Error","attrib4lab2Reducingpricesby15"], r_base_CN$estimates$attrib4lab["Std. Error","attrib4lab1Reducingpricesby30"])
c_attribute6 <- c(0, r_base_CN$estimates$attrib6lab["AMCE","attrib6lab2Occasionalcampaigns"], r_base_CN$estimates$attrib6lab["AMCE","attrib6lab1Frequentcampaigns"])
s_attribute6 <- c(0, r_base_CN$estimates$attrib6lab["Std. Error","attrib6lab2Occasionalcampaigns"], r_base_CN$estimates$attrib6lab["Std. Error","attrib6lab1Frequentcampaigns"])


c <- c(c_attribute1, c_attribute2, c_attribute7, c_attribute3, c_attribute5, c_attribute4, c_attribute6)
s <- c(s_attribute1, s_attribute2, s_attribute7, s_attribute3, s_attribute5, s_attribute4, s_attribute6)

df_r_base_CN <- data.frame(c,s)
df_r_base_CN$att <- c(r_base_CN$attributes$attrib1lab[1], r_base_CN$attributes$attrib1lab[2:3], 
                      r_base_CN$attributes$attrib2lab[1], r_base_CN$attributes$attrib2lab[2:5],
                      r_base_CN$attributes$attrib7lab[3], r_base_CN$attributes$attrib7lab[2], r_base_CN$attributes$attrib7lab[1],
                      r_base_CN$attributes$attrib3lab[5], r_base_CN$attributes$attrib3lab[4], r_base_CN$attributes$attrib3lab[3], r_base_CN$attributes$attrib3lab[2], r_base_CN$attributes$attrib3lab[1], 
                      r_base_CN$attributes$attrib5lab[3], r_base_CN$attributes$attrib5lab[2], r_base_CN$attributes$attrib5lab[1],
                      r_base_CN$attributes$attrib4lab[3], r_base_CN$attributes$attrib4lab[2], r_base_CN$attributes$attrib4lab[1], 
                      r_base_CN$attributes$attrib6lab[3], r_base_CN$attributes$attrib6lab[2], r_base_CN$attributes$attrib6lab[1])
df_r_base_CN$att <- factor(df_r_base_CN$att, levels = c(rev(c("3Nonewtax", "2Increasingpricesby15", "1Increasingpricesby30",
                                                              "2Generalgovernmentbudget", "1Publicenvironmentalandclimateprotectionprograms", "3Publicprogramsforlowincomehouseholds", "4Reduceincometaxes", "5Notaxrevenues",
                                                              "3Keepingsubsidiesatcurrentlevel", "2Halvingsubsidies", "1Eliminatingsubsidies",
                                                              "5Neverrestricted", "41dayperweek", "33daysperweek", "25daysperweek", "1Alwaysrestricted", 
                                                              "3Nonewrequirements", "2Atleast40milespergallon", "1Atleast55milespergallon",
                                                              "3Noincreaseofsupport", "2Reducingpricesby15", "1Reducingpricesby30",
                                                              "3Nocampaigns", "2Occasionalcampaigns", "1Frequentcampaigns"))),
                           labels = c(rev(c("Baseline: No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                            "Baseline: General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                            "Baseline: Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                            "Baseline: Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                            "Baseline: No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                            "Baseline: No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                            "Baseline: No campaigns", "Occasional campaigns", "Frequent campaigns"))))

df_r_base_CN$policy <- ifelse(df_r_base_CN$att == "Baseline: No new tax" |  df_r_base_CN$att == "Increasing prices by 15%" |  df_r_base_CN$att == "Increasing prices by 30%", "1", ifelse(
  df_r_base_CN$att == "Baseline: General government budget" | df_r_base_CN$att == "Public environmental and climate protection programs" | df_r_base_CN$att == "Public programs for low-income households" | df_r_base_CN$att == "Reduce income taxes" | df_r_base_CN$att == "No tax revenues", "2", ifelse(
    df_r_base_CN$att == "Baseline: Keeping subsidies at current level" | df_r_base_CN$att == "Halving subsidies" | df_r_base_CN$att == "Eliminating subsidies", "3", ifelse(
      df_r_base_CN$att == "Baseline: Never restricted" | df_r_base_CN$att == "1 day per week" | df_r_base_CN$att == "3 days per week" | df_r_base_CN$att == "5 days per week" | df_r_base_CN$att == "Always restricted", "4", ifelse(
        df_r_base_CN$att == "Baseline: No new requirements" | df_r_base_CN$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_r_base_CN$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          df_r_base_CN$att == "Baseline: No campaigns" | df_r_base_CN$att == "Occasional campaigns" | df_r_base_CN$att == "Frequent campaigns", "6", ifelse(
            df_r_base_CN$att == "Baseline: No increase of support" | df_r_base_CN$att == "Reducing prices by 15%" | df_r_base_CN$att == "Reducing prices by 30%", "7",
            "8")))))))

df_r_base_CN$policy = factor(df_r_base_CN$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))

p_r_base_CN <- ggplot(df_r_base_CN, aes(x=att, y = c)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s)) +
  theme_bw() + 
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") 
p_r_base_CN


################## US ###############

df_conj_US <- subset(df_conj, country == "US" & incgrp=="4th Quintile" | country == "US" & incgrp=="5th Quintile")
r_base_US <- amce(rate ~ attrib1_lab +  attrib2_lab +  attrib7_lab +  attrib3_lab +  
                    attrib5_lab +  attrib4_lab +  attrib6_lab, data = df_conj_US, 
                  cluster=TRUE, respondent.id="id", 
                  design=design, baselines = baselines_base)


plot(r_base_US)


c_attribute1 <- c(0, r_base_US$estimates$attrib1lab["AMCE","attrib1lab2Increasingpricesby15"], r_base_US$estimates$attrib1lab["AMCE","attrib1lab1Increasingpricesby30"])
s_attribute1 <- c(0, r_base_US$estimates$attrib1lab["Std. Error","attrib1lab2Increasingpricesby15"], r_base_US$estimates$attrib1lab["Std. Error","attrib1lab1Increasingpricesby30"])
c_attribute2 <- c(0, r_base_US$estimates$attrib2lab["AMCE","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], r_base_US$estimates$attrib2lab["AMCE","attrib2lab3Publicprogramsforlowincomehouseholds"],r_base_US$estimates$attrib2lab["AMCE","attrib2lab4Reduceincometaxes"], r_base_US$estimates$attrib2lab["AMCE","attrib2lab5Notaxrevenues"])
s_attribute2 <- c(0, r_base_US$estimates$attrib2lab["Std. Error","attrib2lab1Publicenvironmentalandclimateprotectionprograms"], r_base_US$estimates$attrib2lab["Std. Error","attrib2lab3Publicprogramsforlowincomehouseholds"],r_base_US$estimates$attrib2lab["Std. Error","attrib2lab4Reduceincometaxes"],r_base_US$estimates$attrib2lab["Std. Error","attrib2lab5Notaxrevenues"])
c_attribute7 <- c(0, r_base_US$estimates$attrib7lab["AMCE","attrib7lab2Halvingsubsidies"], r_base_US$estimates$attrib7lab["AMCE","attrib7lab1Eliminatingsubsidies"])
s_attribute7 <- c(0, r_base_US$estimates$attrib7lab["Std. Error","attrib7lab2Halvingsubsidies"], r_base_US$estimates$attrib7lab["Std. Error","attrib7lab1Eliminatingsubsidies"])
c_attribute3 <- c(0, r_base_US$estimates$attrib3lab["AMCE","attrib3lab41dayperweek"], r_base_US$estimates$attrib3lab["AMCE","attrib3lab33daysperweek"], r_base_US$estimates$attrib3lab["AMCE","attrib3lab25daysperweek"], r_base_US$estimates$attrib3lab["AMCE","attrib3lab1Alwaysrestricted"])
s_attribute3 <- c(0, r_base_US$estimates$attrib3lab["Std. Error","attrib3lab41dayperweek"], r_base_US$estimates$attrib3lab["Std. Error","attrib3lab33daysperweek"], r_base_US$estimates$attrib3lab["Std. Error","attrib3lab25daysperweek"], r_base_US$estimates$attrib3lab["Std. Error","attrib3lab1Alwaysrestricted"])
c_attribute5 <- c(0, r_base_US$estimates$attrib5lab["AMCE","attrib5lab2Atleast40milespergallon"], r_base_US$estimates$attrib5lab["AMCE","attrib5lab1Atleast55milespergallon"])
s_attribute5 <- c(0, r_base_US$estimates$attrib5lab["Std. Error","attrib5lab2Atleast40milespergallon"], r_base_US$estimates$attrib5lab["Std. Error","attrib5lab1Atleast55milespergallon"])
c_attribute4 <- c(0, r_base_US$estimates$attrib4lab["AMCE","attrib4lab2Reducingpricesby15"], r_base_US$estimates$attrib4lab["AMCE","attrib4lab1Reducingpricesby30"])
s_attribute4 <- c(0, r_base_US$estimates$attrib4lab["Std. Error","attrib4lab2Reducingpricesby15"], r_base_US$estimates$attrib4lab["Std. Error","attrib4lab1Reducingpricesby30"])
c_attribute6 <- c(0, r_base_US$estimates$attrib6lab["AMCE","attrib6lab2Occasionalcampaigns"], r_base_US$estimates$attrib6lab["AMCE","attrib6lab1Frequentcampaigns"])
s_attribute6 <- c(0, r_base_US$estimates$attrib6lab["Std. Error","attrib6lab2Occasionalcampaigns"], r_base_US$estimates$attrib6lab["Std. Error","attrib6lab1Frequentcampaigns"])


c <- c(c_attribute1, c_attribute2, c_attribute7, c_attribute3, c_attribute5, c_attribute4, c_attribute6)
s <- c(s_attribute1, s_attribute2, s_attribute7, s_attribute3, s_attribute5, s_attribute4, s_attribute6)

df_r_base_US <- data.frame(c,s)
df_r_base_US$att <- c(r_base_US$attributes$attrib1lab[1], r_base_US$attributes$attrib1lab[2:3], 
                      r_base_US$attributes$attrib2lab[1], r_base_US$attributes$attrib2lab[2:5],
                      r_base_US$attributes$attrib7lab[3], r_base_US$attributes$attrib7lab[2], r_base_US$attributes$attrib7lab[1],
                      r_base_US$attributes$attrib3lab[5], r_base_US$attributes$attrib3lab[4], r_base_US$attributes$attrib3lab[3], r_base_US$attributes$attrib3lab[2], r_base_US$attributes$attrib3lab[1], 
                      r_base_US$attributes$attrib5lab[3], r_base_US$attributes$attrib5lab[2], r_base_US$attributes$attrib5lab[1],
                      r_base_US$attributes$attrib4lab[3], r_base_US$attributes$attrib4lab[2], r_base_US$attributes$attrib4lab[1], 
                      r_base_US$attributes$attrib6lab[3], r_base_US$attributes$attrib6lab[2], r_base_US$attributes$attrib6lab[1])
df_r_base_US$att <- factor(df_r_base_US$att, levels = c(rev(c("3Nonewtax", "2Increasingpricesby15", "1Increasingpricesby30",
                                                              "2Generalgovernmentbudget", "1Publicenvironmentalandclimateprotectionprograms", "3Publicprogramsforlowincomehouseholds", "4Reduceincometaxes", "5Notaxrevenues",
                                                              "3Keepingsubsidiesatcurrentlevel", "2Halvingsubsidies", "1Eliminatingsubsidies",
                                                              "5Neverrestricted", "41dayperweek", "33daysperweek", "25daysperweek", "1Alwaysrestricted", 
                                                              "3Nonewrequirements", "2Atleast40milespergallon", "1Atleast55milespergallon",
                                                              "3Noincreaseofsupport", "2Reducingpricesby15", "1Reducingpricesby30",
                                                              "3Nocampaigns", "2Occasionalcampaigns", "1Frequentcampaigns"))),
                           labels = c(rev(c("Baseline: No new tax", "Increasing prices by 15%", "Increasing prices by 30%", 
                                            "Baseline: General government budget", "Public environmental and climate protection programs", "Public programs for low-income households", "Reduce income taxes", "No tax revenues",
                                            "Baseline: Keeping subsidies at current level", "Halving subsidies", "Eliminating subsidies",
                                            "Baseline: Never restricted", "1 day per week", "3 days per week", "5 days per week", "Always restricted", 
                                            "Baseline: No new requirements", "Maximum consumption of 5.9 litres per 100 kilometres", "Maximum consumption of 4.3 litres per 100 kilometres",
                                            "Baseline: No increase of support", "Reducing prices by 15%", "Reducing prices by 30%",
                                            "Baseline: No campaigns", "Occasional campaigns", "Frequent campaigns"))))

df_r_base_US$policy <- ifelse(df_r_base_US$att == "Baseline: No new tax" |  df_r_base_US$att == "Increasing prices by 15%" |  df_r_base_US$att == "Increasing prices by 30%", "1", ifelse(
  df_r_base_US$att == "Baseline: General government budget" | df_r_base_US$att == "Public environmental and climate protection programs" | df_r_base_US$att == "Public programs for low-income households" | df_r_base_US$att == "Reduce income taxes" | df_r_base_US$att == "No tax revenues", "2", ifelse(
    df_r_base_US$att == "Baseline: Keeping subsidies at current level" | df_r_base_US$att == "Halving subsidies" | df_r_base_US$att == "Eliminating subsidies", "3", ifelse(
      df_r_base_US$att == "Baseline: Never restricted" | df_r_base_US$att == "1 day per week" | df_r_base_US$att == "3 days per week" | df_r_base_US$att == "5 days per week" | df_r_base_US$att == "Always restricted", "4", ifelse(
        df_r_base_US$att == "Baseline: No new requirements" | df_r_base_US$att == "Maximum consumption of 5.9 litres per 100 kilometres" | df_r_base_US$att == "Maximum consumption of 4.3 litres per 100 kilometres", "5", ifelse(
          df_r_base_US$att == "Baseline: No campaigns" | df_r_base_US$att == "Occasional campaigns" | df_r_base_US$att == "Frequent campaigns", "6", ifelse(
            df_r_base_US$att == "Baseline: No increase of support" | df_r_base_US$att == "Reducing prices by 15%" | df_r_base_US$att == "Reducing prices by 30%", "7",
            "8")))))))

df_r_base_US$policy = factor(df_r_base_US$policy, levels=c(1:8), labels=c("Tax","Earmarking","Industry Subsidy", "Access restriction", "Standards", "Campaign", "PT prices", "Other"))

p_r_base_US <- ggplot(df_r_base_US, aes(x=att, y = c)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s)) +
  theme_bw() + 
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") 
p_r_base_US


####### Rating overview (with 95% confidence interval)

#country
summary(df_conj$rate[df_conj$country=="CN" & df_conj$incgrp=="4th Quintile" | df_conj$country=="CN" & df_conj$incgrp=="5th Quintile"])
summary(df_conj$rate[df_conj$country=="DE" & df_conj$incgrp=="4th Quintile" | df_conj$country=="DE" & df_conj$incgrp=="5th Quintile"])
summary(df_conj$rate[df_conj$country=="US" & df_conj$incgrp=="4th Quintile" | df_conj$country=="US" & df_conj$incgrp=="5th Quintile"])
rate_summary_country <- summarySE(df_conj, measurevar="rate", groupvars=c("country"), na.rm = T, conf.interval = 0.95)
rate_summary_country

fig_rate_summary_country1 <- ggplot(rate_summary_country, aes(x=country, y=rate, fill=country)) + 
  geom_bar(position=position_dodge(), stat="identity",
           colour="black",
           size=0.3) + 
  geom_errorbar(aes(ymin=rate-ci, ymax=rate+ci),
                size=.3,
                width=.2,
                position=position_dodge(.9)) +
  theme_bw() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        panel.grid.major.y = element_line(size=.5, color="gray80" )) +
  xlab(" ") +
  ylab("Average rating") +
  #theme(plot.title = element_text(size = 20, face = "bold"),axis.text.x =element_text(size=20, face ="bold"),axis.text.y =element_text(size=20, face ="bold"),strip.text.x = element_text(size=20, face ="bold"), axis.title.y = element_text(size = 20, face ="bold"))+
  #ggtitle("Average Rating (1-7) by Country") +
  scale_fill_manual("", values = c("US" = "#b2abd2", "DE" = "#5e3c99", "CN" = "#e66101"), labels=c("US" = "USA", "DE" = "Germany", "CN" = "China")) +
  scale_y_continuous(breaks = c(1,2,3,4,5,6,7), labels = c("1","2", "3", "4", "5", "6", "7")) +
  scale_x_discrete(labels = c("China", "Germany", "USA")) +
  coord_cartesian(ylim = c(1.27, 6.73)) 

legend3 <- get_legend(fig_rate_summary_country1)

fig_rate_summary_country <- fig_rate_summary_country1 + guides(fill=FALSE)

#overview average rating
df_conj_45 <- subset(df_conj, df_conj$incgrp=="4th Quintile" | df_conj$incgrp=="5th Quintile")

rateavrg <- aggregate(rate ~ attrib1_lab +  attrib2_lab +  attrib3_lab +  attrib4_lab +  
            attrib5_lab +  attrib6_lab +  attrib7_lab + country, df_conj_45, mean)

fig_rate_average <- ggplot(rateavrg, aes(x=rate, colour=country)) + 
  geom_density(size=1.2) + 
  theme_bw() +
  theme(panel.grid.major.x = element_line(size=.5, color="gray80"),
        panel.grid.minor.y = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.5, color="gray80" )) +
  xlab("Distribution of average ratings for policy-packages") +
  ylab("Density") +
  scale_colour_manual(values = c("US" = "#b2abd2", "DE" = "#5e3c99", "CN" = "#e66101")) +
  scale_x_continuous(breaks = c(1,2,3,4,5,6,7), labels = c("1","2", "3", "4", "5", "6", "7")) +
  coord_cartesian(ylim = c(.0257, 0.55), xlim = c(1.27, 6.73)) +
  #geom_vline(xintercept = 4.740764, color = "#e66101", size=1.2, linetype = "longdash") +
  #geom_vline(xintercept = 3.651830, color = "#5e3c99", size=1.2, linetype = "longdash") +
  #geom_vline(xintercept = 4.069573, color = "#b2abd2", size=1.2, linetype = "longdash") +
  guides(colour=FALSE)



#############paper graphs export


#conjoint graph by country, highest income classes
df_countries <- rbind(df_r_base_US, df_r_base_CN, df_r_base_DE, df_c_base_US, df_c_base_CN, df_c_base_DE)
df_countries$ctr <- rep(c("USA", "China", "Germany", "USA", "China", "Germany"), each=nrow(df_r_base_US))
df_countries$grp <- rep(c("Rating", "Rating", "Rating", "Choice", "Choice", "Choice"), each=nrow(df_r_base_US))


cntrs_output3 <- ggplot(df_countries, aes(x=att, y = c, group = ctr, shape = ctr, colour = ctr)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s), position = position_dodge(width=.5)) +
  theme_bw() + 
  theme(panel.grid.major.x = element_line(size=.3, color="gray80"),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_line(size=.3, color="gray80" ),
        panel.grid.minor.y = element_blank()) +
  coord_flip() +
  ylab("Estimated AMCE") +
  xlab("Attributes") +
  #ggtitle("Estimated Average Marginal Component Effects") +
  geom_hline(yintercept = 0, lty="dashed") + 
  facet_grid(policy~grp, scales = "free") +
  scale_color_manual(name="Country", values=c("#e66101", "#5e3c99", "#b2abd2")) +
  scale_shape_discrete(name="Country") + 
  theme(strip.background = element_rect(fill = 'white'))

cntrs_output2 <- cntrs_output3 + theme(legend.position="none")

legend <- get_legend(cntrs_output3)


graph21 <- ggdraw() +
  draw_plot(fig_rate_summary_country, x = 0, y = 0, width = .25, height = .25) +
  draw_plot(fig_rate_average, x = .25, y = 0, width = .6, height = .25) +
  draw_plot(cntrs_output2, x = 0, y = 0.25, width = 1, height = 0.75) +
  draw_plot(legend3, x = 0.85, y = 0.065, width = 0.15, height = 0.125) +
  draw_plot(legend, x = 0.85, y = 0.14, width = 0.15, height = 0.125) +
  draw_plot_label(label = c("A", "B", "C"), size = 14,
                  x = c(0, 0, 0.25), y = c(1, 0.25, 0.25))
graph21

ggsave(filename = "./figs/graph21_quintile45.png", graph21, width = 9, height = 11)