b-lux/RTutorCarbonLeakage
RTutor problem set Optimized ndustry compensation under the threat of carbon leakage

Package index
Search the b-lux/RTutorCarbonLeakage package
Vignettes
Functions
0
Source code
2
Man pages
0
Home
/
GitHub
/
b-lux/RTutorCarbonLeakage
/

In b-lux/RTutorCarbonLeakage: RTutor problem set Optimized ndustry compensation under the threat of carbon leakage

Problemset Optimized industry compensation under the threat of carbon leakage

Author: Benjamin Lux

< ignore

Run code below to generate and show the problemset

library(RTutor)
library(yaml)

setwd("/Users/Benny/Dropbox/Documents/Uni/Abschlussarbeiten/Masterarbeit/R-Tutor/RTutorCarbonLeakage") #adapt path
ps.name = "Optimal_permit_allocation_carbon_leakage"
sol.file = paste0(ps.name, "_sol.Rmd")
libs = c("ggplot2", "dplyr", "tidyr", "foreign", "pracma", "regtools", "yarrr")

# Generate problem set files in working directory
create.ps(sol.file = sol.file, ps.name = ps.name, user.name = NULL,libs = libs, extra.code.file = "fun.r", var.txt.file = "Optimal_permit_allocation_carbon_leakage_Var.txt", stop.when.finished = FALSE, addons="quiz")
# 
# Show in webbrowser. You can adapt the argumets below
show.ps(ps.name, load.sav = FALSE, launch.browser = TRUE, sample.solution = FALSE, is.solved = FALSE)

>

In their paper "Industry Compensation under Relocation Risk: A Firm-Level Analysis of the EU Emissions Trading Scheme" Ralf Martin, Mirabelle Muûls, Laure B. de Preux and Ulrich J. Wagner (2014) analyzed how to distribute free of charge $CO_2$ emission certificates optimally among firms to encounter the threat of carbon leakage. In the following, I will refer to this paper by Martin et al. The aim of this interactive R tutorial is to reproduce and discuss their findings.
The website of the American Economic Association provides the original paper and corresponding data. This RTutor problem can be accessed here.

Exercise Overview

The main interest of Martin et al.'s paper was to find an optimal solution to a problem called carbon leakage that emerges as $CO_2$ emissions are imposed with a price and an overall cap under the European Union Emissions Trading System (EU ETS). Firms might shift their greenhouse gas (GHG) emissions to unregulated places outside the European Union (EU) in order to avoid associated costs. This phenomenon is referred to as carbon leakage (Commission Decision 2010/2/EU). There are two main problems with this response to putting a price tag on carbon emissions. First of all, it weakens the effectiveness of the environmental regulation. Since $CO_2$ is a global public bad, there is no advantage in a regulation trying to reduce overall contribution to global warming, but in reality only shifting emissions somewhere else. Second, Martin et al. argued that it might have a bad economic impact on the EU as well, since relocating emissions goes hand in hand with relocating jobs and taxable profits. Policy measures in place offer compensation in form of free permits to firms that are either carbon intensive and/or trade exposed in order to prevent them from relocation. Martin et al. showed that this approach of compensating the firms is less efficient than their own proposal of an alternative compensation rule. They found that the aggregate amount of jobs at risk to be relocated can be reduced by more than half without extending the compensation level. In the course of this optimization the authors used data on firms' vulnerability to carbon pricing that usually is not available to policy-makers. Therefore, they based their compensation rule in a next step on publicly observable firm data. Even with these easy compensation rules the risk of job relocation can be substantially reduced.

Structure of the problemset:

The problemset provides an interactive environment to reproduce and discuss Martin et al.'s core findings and is structured in exercises. The interaction evolves as you, the user, have the possibility to write little R code fragments to establish a basis for economic conclusions. Guidance is provided along the way. Answer the short quizzes to deepen your understanding of the content and to test your economic intuition. It is not obligatory, but recommended, to solve the exercises in the given order. R-skills acquired in the first exercises are assumed to be known in the further course of the problemset. If you are primarily interested in the economic content of the paper you might want to skip Exercise 3.2, which gives a technical, yet illustrative, introduction into dynamic programming.

Exercise 1 Current Allocation Scheme

We start our journey into the EU ETS by introducing some key terms that are necessary to discuss its features. Afterwards we will visualize the prevailing policy measures aiming at the prevention of carbon leakage and attempt a first evaluation.

i) The initial idea

It all started with the insight that climate change is really happening and that anthropogenic GHG emissions are a driving factor. Today, there is no serious scientific doubt that human activities such as burning fossil fuels in transport and industrial processes, chopping down forests and animal husbandry increase the concentration of greenhouse gases in the atmosphere and thus are responsible for global warming (Stern, 2007). In his review Stern (2007) claimed that the cost of the consequences of global warming, including floods and droughts, will exceed the cost of early prevention measures. To limit the consequences of global warming legally-binding GHG reduction targets for 37 industrialized countries were established in the Kyoto Protocol in 1997. Within the Kyoto framework the emission reduction target of the EU was set to 8% compared to the level in 1990 in the first commitment period (2008 - 2012). Subsequently it became necessary to develop union-wide policy instruments that could guarantee the committed reduction of emissions. The cornerstone of the EU's reduction approach is the EU ETS (European Commission, 2016). As a cap and trade system it sets a limit to GHG emissions from stationary sources - mostly power stations and industrial plants - and aircraft operators. The overall cap ensures compliance with Kyoto emission reduction targets. This cap is translated into pollution permits, which participants have to acquire to emit $CO_2$. Permits can be traded amongst emitters. The trading approach results in a carbon price that shall incentivize participants to cut emissions in the most cost-effective manner (European Commission, 2016).

Visualization of the different EU ETS initialization stages. Figure 1.1: EU ETS initialization stages. Source: Author's own graph.

As visualized in Figure 1.1, the installation of the EU ETS was realized in different stages. Phases I and II lasting from 2005 until 2012 can be seen as test phases, since almost all permits were allocated to firms for free. The amount of emission allowances allocated to firms was "based on historical emissions and adjusted for growth projections and the national contribution towards the EU's joint emission target" (Martin et al, 2014, p. 2486). We will refer to this allocation method as grandfathering. More detailed information on grandfathering are given in the info-box below.

< info "Grandfathering"

The grandfathering approach is characterized in the EU ETS Handbook (European Commission, 2016) as follows:

Even though the member states were allowed to auction off up to 5% of allowances in phase I and up to 10% in phase II, Ellerman and Joskow (2008) found that auctioning was utilized rarely.

>

Trading phase III, lasting from 2013 until 2020, is characterized by the gradual transition from grandfathering towards full auctioning as the default allocation principle (European Commission, 2016). The auctioning of permits fits the polluter-pays principle, a corner stone of environmental policy-making in the EU (Directive 2004/35/CE). Martin et al. noted that the principle of auctioning constitutes transfer payments from polluters to governments and taxpayers. The transition is realized by a benchmarking-approach (Commission Decision 2011/87/EU). The amount of free permits allocated to firms is now based on product-related benchmarks and decreases gradually over time. As firms receive less pollution allowances for free step by step, they need to acquire missing permits in auctions or reduce their carbon emissions altogether. It is envisioned that by 2020 the share of free permits is reduced to only 30% of the respective benchmark (European Commission, 2016). See the info-box below for more detailed information on benchmarking.

Note: The EU ETS has constantly increased in terms of geography, sectors and GHGs since its beginning in 2005. In the course of this problem set we are going to concentrate on manufacturing industries and $CO_2$ emissions, since these are both the most important and most interesting topics.

< info "Benchmarking"

Under benchmarking the amount of free permits $q_{ijt}^b$ an installation $i$ manufacturing a product $j$ receives in year $t$ is determined by the following formula (Martin et al., 2014, p. 2486):

$$ q_{ijt}^b = benchmark_j \cdot historical\: activity\: level_{i,j} \cdot reduction_{j,t} \cdot correction_t. $$

>

< quiz "Introduction of the EU ETS"

parts: - question: 1. What was the EU's GHG reduction target (in % compared to 1990) for the first commitment period of the Kyoto Protocol?
answer: 8 roundto: 1

>

If you want to find out more about the development of the EU ETS have a look at the ETS Summer University. The European Commission has set up a comprehensive online course on the background and design of Emission Trading Systems, here.

ii) The Carbon Leakage Decision

It comes with no surprise that carbon intensive industries complained that the increased cost resulting from full auctioning of emission permits were a competitive disadvantage compared to competitors in less regulated jurisdictions. This implies the threat that companies could transfer their production to countries imposing laxer emission standards to save carbon costs. Since carbon emissions are an "international externality" (Markusen, 1975) the European Commission (EC) recognized that this policy "could lead to an increase in greenhouse gas emissions in third countries where industry would not be subject to comparable carbon constraints (‘carbon leakage’) and undermine the environmental integrity and benefit of actions by the Union” (Commission Decision 2010/2/EU, p.1). Moreover, it costs jobs and taxable profits in the EU. The EC gave in to this threat and guaranteed 100% of benchmark allocations at no cost to sectors deemed to be at risk of carbon leakage in the Carbon Leakage Decision (Commission Decision 2010/2/EU). According to this Commission Decision carbon intensity (CI) and trade intensity (TI) are the relevant quantities that determine a sector's leakage risk. Have a look at the next info-box for a formal definition of these terms. Based on these measures the EC defined thresholds to determine sectors or sub-sectors to be exempt from permit auctioning.

< info "Carbon and Trade Intensity"

$$ CI = \frac{( direct\: CO_2\: emissions + indirect\: CO_2\: emissions)\cdot CO_2\: price}{gross\: value\: added\: of\: a\: sector} $$

CI represents the additional production costs imposed by compliance with the EU ETS. It is calculated as the sum of the direct and indirect cost of full permit auctioning as a proportion of gross value added. The direct costs are determined by multiplying direct $CO_2$ emissions with a carbon price of €30/tCO2. The indirect costs account for higher electricity prices, that occur since the power generation sector is not exempt from full auctioning. These indirect costs are the product of electricity consumption (in MWh), the average $CO_2$ emissions of power generation in the EU (0.465 tCO2/MWh), and the carbon price (€30/tCO2).

$$ TI = \frac{non\:EU\: Exports +non\:EU\: Imports}{Market\: size} = \frac{non\:EU\: Exports +non\:EU\: Imports}{Annual\: turnover + non\:EU\: Imports} $$ TI is "defined as the ratio between the total value of exports to third countries plus the value of imports from third countries and the total market size for the Union (annual turnover plus total imports from third countries)" (Commission Decision 2010/2/EU, p.1).

>

In the following, we are using firm-level data to visualize the exemption criteria established under the Carbon Leakage Decision. This will help to understand and to analyze the policy measures in place to reduce the relocation risk.

All the data frames being used throughout this problem set are contained in a folder called data. In order to work with a data frame it needs to be loaded into the workspace of our problem set.

< info "read.table()"

Tabulated data stored in a file called dataframe can be loaded to the workspace using the read.table() command. The file name of the data frame of interest is expressed in relative terms of the current working directory and needs to be enclosed in quotation marks. As in the example below, this dataframe can be assigned to a new variable data.

data = read.table("data/dataframe")

>

Task: In this example the read.table() command is used to load the data frame firm.data and assign it to a variable called firm.data. The second command, head(), prints the first entries of the data frame. Execute the code by hitting the check-button. Further information on the read.table() command can be found in the info-box above.

#< task
firm.data = read.table("data/firm.data")
head(firm.data)
#>

This data set was prepared beforehand for the purpose of this exercise to start quickly with some interesting graphs. We will get to know the tools being used in the preparation process in the course of this and following exercises. The interested reader is referred to the footnote of this exercise to see the data preparation code. In the data set each row displays information on one firm, e.g. to which sector it belongs, its carbon or trade intensity or its number of employees. To have a look at the whole data set click on the data-button above the command window. More detailed information on the variables and the data frame are given in the info-box below. Short variable descriptions are available by hovering over the variable names with your cursor, too.

< info "firm.data"

The data frame contains information on 3814 manufacturing firms. The figures on carbon and trade intensity were calculated on the data basis of the Community Independent Transactions Log (CITL) and EUROSTAT. The data on $CO_2$ emissions stem from CITL too. Employment figures are obtained from the ORBIS database maintained by Bureau Van Dijk.

>

< info "NACE"

NACE is an acronym for Nomenclature statistique des activités économiques dans la Communauté européenne. It provides a systematic classification of productive economic activities. Its purpose is to improve comparability of economic statistics like employment or production. The hierarchical structure of NACE classifies economic sectors up to a four digit code. The code 1752 for example refers to manufacturing of cordage, rope, twine and netting. The full mapping of the codes is provided here.

>

First, we will take an overview of how sectors (NACE four-digit) relate to the EC's definitions of carbon and trade intensity. Therefore, we need to transform the firm-level data into sectoral data. When it comes to data manipulation the dplyr package is a good choice. It is powerful and intuitive to read and write. To aggregate information on the sector level we will use a combination of two commands that we will encounter quite frequently: group_by() and summarize().

< info "group_by(), summarize() and the 'pipe' "

The dplyr package written by Wickham and Francois (2016) provides nice tools to work with data frames.

dataframe = group_by(dataframe, var)

summarize(dataframe, 
          mean_1 = mean(var1),
          median_2 = median(var2))

dataframe %>% group_by(var2) %>% summarize(mean(var1))

A very useful cheat sheet on the dplyr package is available here.

>

Task: Delete the #-signs and complete the code in the chunk below by replacing ???. Calculate summarized figures on the sector level (sec4dig) like the average carbon intensity or the aggregate number of employees. Examples on how to use groupe_by() and summarize() are provided in the info-box above. Results are stored in a new data frame called fig1.2.data. A hint on how to solve this task can be obtained by clicking the hint-button. Hit check to execute the code once you are done. If you cannot or do not want to solve the task yourself click solution and check afterwards.

#< task_notest
# Load the required package 'dplyr' from the library:
library(dplyr)

# Complete the code below by replacing ??? with your solution

# fig1.2.data = firm.data %>%
#   group_by( ??? ) %>%
#   summarize(carb_int = mean(carb_int, na.rm = TRUE),
#             trade_int = mean(trade_int, na.rm = TRUE),
#             n_install = sum(n_install, na.rm = TRUE),
#             emp = sum(emp, na.rm = TRUE),
#             co2 = sum(co2, na.rm = TRUE),
#             firm_sec = sum(firm_sec, na.rm=TRUE))
#>

fig1.2.data = firm.data %>%
  group_by(sec4dig) %>%
  summarize(carb_int = mean(carb_int, na.rm = TRUE),
            trade_int = mean(trade_int, na.rm = TRUE),
            n_install = sum(n_install, na.rm = TRUE),
            emp = sum(emp, na.rm = TRUE),
            co2 = sum(co2, na.rm = TRUE),
            firm_sec = sum(firm_sec, na.rm=TRUE))
#< hint
display("Replace ??? with 'sec4dig'.")
#>

These aggregate information are now used to generate the first graph on our own. For this purpose, we are going to use the package ggplot2, which provides us with a great variety of tools to create nice plots.

< info "ggplot2"

After dplyr, ggplot2 is the second package written by Wickham (2009) that we encounter. Since it enables us to compose a graph by combining many independent components, it is very powerful. Graph creation takes place in multiple layers (Wickham, 2009). In the first layer we can specify the data dataframe (containing variables var1 and var2) and default x- and y-variables we want to work with.

fig = ggplot(data = dataframe, aes(x = var1, y = var2))

In the next layer we choose the type of graph. Different layers are combined with the +-operator. If we want to display our data with points we add geom_point().

fig = fig + geom_point()

In this manner many different layers can be added to the graph, including regression lines or error bars.

Many examples on plot creation with ggplot2 can be found here.

>

Task: Create a bubble graph in which the horizontal axis shows a sector's trade intensity trade_int and the vertical axis its carbon intensity carb_int. A third dimension is added by setting the size of the data points equal to the number of installations n_install in a sector. The required data is stored in the previously summarized data frame fig1.2.data. Parts of the code are already provided. Help on how to use ggplot2 can be found in the info-section above. As before, further help is provided by clicking the hint-button. Complete the code and press the check-button to create and display the graph. 

#< task_notest
# Load the required package 'ggplot2' from the library:
library(ggplot2)

# The argument 'shape=1' in 'geom_point' displays data 
# points as hollow circles. 'scale_size_area()' ensures 
# that we scale area instead of radius for size.

# fig1.2 = ggplot(data = ??? , aes(x = ??? , y = ??? )) +
#        geom_point(aes(size = ??? ), shape=1) +
#        scale_size_area(max_size=22) +
#        xlab("Trade Intensity") +
#        ylab("Carbon Intensity") 

# Show the plot:
# fig1.2
#>

fig1.2 = ggplot(data = fig1.2.data, aes(x=trade_int, y=carb_int)) +
       geom_point(aes(size = n_install), shape=1) +
       scale_size_area(max_size=22) +
       xlab("Trade Intensity") +
       ylab("Carbon Intensity") 
fig1.2
#< hint
display("Do not forget to remove the #-signs in front of the code you want to execute. The assignments you have to do are 'data = fig1.2.data', 'x = trade_int', 'y=carb_int' and 'size=n_install'.")
#>

Figure 1.2: Distribution of manufacturing industries with respect to CI and TI. Source: Adapted from Martin et al. (2014, p. 2488).

Note: Martin et al.'s graph in the paper differs slightly from this one because they did not exclude non-ETS installations from their sample. Since the graph is supposed to give insight into the EU ETS and its exemption criteria, I altered the data preparation process in this aspect.

Figure 1.2 illustrates how the different manufacturing industries are distributed with respect to CI (vertical axis) and TI (horizontal axis). Each circle represents one industry and its size (area) is proportional to the number of installations contained. We can observe a high density of small sectors, containing less than 200 installations, with a CI below 5% and a TI between 0% and 90%. The most medium- to large-sized sectors with respect to number of installations contained, display a CI between 5% and 20% and a TI between 0% and 40%. Only four sectors have a CI greater than 40%. None of the sectors displays a very high CI in combination with a very high TI.

It comes handy that new geometric objects can be added to an existing graph easily with the + -operator under ggplot2. This helps us to visualize the EC's definition of "sectors or sub sectors deemed to be exposed to a significant risk of carbon leakage" (Commission Decision 2010/2/EU, p.1). According to the EC significant relocation risk of a sector is given if the following thresholds for carbon or trade intensity are exceeded: $CI > 5\,\%\: \&\: TI > 10\,\%$, or $CI > 30\,\%$, or * $TI > 30\,\%$.

Firms in sectors exceeding these thresholds will continue to receive 100% of benchmark allocations for free. Martin et al. identified three mutually exclusive groups of sectors meeting the stipulated requirements. These are:

Each of these groups can be visualized by a rectangle in our graph.

Task: Following the example for category $A$ given below, fill in the respective values for categories $B$ and $C$ to visualize them as rectangles. Enter a command that shows the modified graph afterwards. The hint-button provides help on how to solve the problem if needed.

#< task_notest 
# Any number of geoms can be added to a ggplot object using the '+' sign. 
# Rectangles can be created by calling 'geom_rect()' and 
# specifying the locations of the four corners.
# 'fill' changes the interior colouring of geometric object.
# 'alpha' sets the transparancy of an object.
# 'annotate()' adds the category names as text.

# fig1.3 = fig1.2 +
#   geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 30, ymax = Inf),   # A
#             fill = "#FBE814", alpha = 0.009) +
#   geom_rect(aes(xmin = ???, xmax = ???, ymin = ???, ymax = ???),   # B
#             fill = "#0072B2", alpha = 0.005) +
#   geom_rect(aes(xmin = ???, xmax = ???, ymin = ???, ymax = ???),   # C
#             fill = "#D73200", alpha = 0.009) +
#   annotate("text", x=c(95, 95,20), y=c(40,25,25), 
#            label = c("A", "B", "C"), size = 10)

#> 

fig1.3 = fig1.2 +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = 30, ymax = Inf), 
            fill = "#FBE814", alpha = 0.009) +
  geom_rect(aes(xmin = 30, xmax = Inf, ymin = -Inf, ymax = 30),
            fill = "#0072B2", alpha = 0.005) +
  geom_rect(aes(xmin = 10, xmax = 30, ymin = 5, ymax = 30),
            fill = "#D73200", alpha = 0.009) +
  annotate("text", x=c(95, 95,20), y=c(40,25,25), 
           label = c("A", "B", "C"), size = 10)

fig1.3

#< hint
display("Take a look at the given example. Following the syntax of the given code fill in the values of the four corners of the rectangles. Add 'fig1' in the end to display the modified graph.")
#>

Figure 1.3: Exempted sectors according to the EC. Source: Adapted from Martin et al. (2014, p. 2488).

Have a closer look on Figure 1.3 and answer the following question:

< quiz "Carbon Leakage Decision"

question: Which exemption category contains the most sectors? sc: - A (yellow) - B (blue)* - C (red) success: Great, your answer is correct! failure: Try again.

>

We can make three interesting observations in Figure 1.3: First, the number of exempt sectors and installations constitutes a relatively large share of the sample. Second, the comparison of the number of circles in the different rectangles reveals that most of the sectors exempted from full auctioning by the EC belong to category $B$. Third, most of these category $B$ sectors have a carbon intensity well below 5% and can therefore not be considered as carbon intensive.

These observations incentivized Martin et al. to dig deeper and have a closer look on the different exempted groups. Following their observations, they split exemption category $B$ into low (CI < 5%) and moderate (5% < CI < 30%) carbon intensity sub-categories. In a next step, we add a new variable to our data frame specifying the exemption category each sector belongs to. The dplyr package provides the required tools: mutate() and filter().

< info "mutate() and filter()"

mutate(data, var2 = 2*var1)

filter(dataframe, var1 > 1)

>

Task: Hit the check-button to add the new variable exempt_group to fig1.2.data. The result is stored in exemption.data and will only contain sectors that could be assigned to one of the exemption categories (use of filter()).

#< task
# The new variable 'exempt_group' is added to the frame. Its values depend on 
# the values of 'carb_int' and 'trade_int' and are selected by
# ifelse(logical condition, TRUE, FALSE).
exemption.data = fig1.2.data %>%
  mutate(exempt_group = 
    ifelse(carb_int>30, "A",
        ifelse(trade_int>30, ifelse(carb_int>5, "B & CI > 5", "B & CI < 5"),
            ifelse(carb_int>5 & trade_int>10, "C","Not exempt")))) %>%
  filter(!is.na(exempt_group))
#>

To learn more about the structure of the different exemption categories, we want to know how firms, employees and $CO_2$ emissions are shared out between them. Therefore, we calculate the share each sector takes in terms of firms, employees and $CO_2$ emissions. Instead of adding new variables to the data frame, we only want to keep the new shares and the exemption category. This can be done with transmute() from the dplyr package. The syntactic structure to call transmute() is equal to the one used for mutate().

Task: Click check to create a data frame share.data which stores the exemption categories and the share of firms, employees and emissions each sector takes in these categories.

#< task_notest

# From the old data frame we only keep the variable 'exempt_group'.
# Besides we only have the share of each sector in the different 
# quantity of interest.
share.data = exemption.data %>%
  transmute(exempt_group,
            share_firms = firm_sec/sum(firm_sec, na.rm = TRUE),
            share_employees = emp/sum(emp, na.rm = TRUE),
            share_emissions = co2/sum(co2, na.rm = TRUE))
#>

Since we want to compare the different exemption categories, we again need to group our data by the exemption groups and sum up the newly introduced shares.

Task: Following the example given above, fill in the missing parts of the code. 

#< task_notest
# The 'summarize_each()' command works in principle similar to the 
# 'summarize()' command. The function called as an argument of 'funs()'
# is applied to all columns but the grouping variable (exempt_group).

# share.data = share.data %>%
#   group_by( ??? ) %>%
#   summarize_each(funs( ??? ))
# share.data

#>
share.data = share.data %>%
  group_by(exempt_group) %>%
  summarize_each(funs(sum))
share.data

#< hint
display("We want to group by 'exempt_group', and 'sum' all the shares in a group up.")
#>

Before we analyze these results, let us display them in a graph. In order to generate a bar graph from this data with ggplot2, we need to transform it from wide format to long format. This means that we gather column names into rows. The tidyr package provides the required function gather().

< info "gather()"

The gather() command takes a data frame dataframe and stores the old column names var1 and var2 under the new key-variable var3. The values previously stored under the old column names are stored in the column value. The source columns that shall be transformed are added separated by comma.

gather(data = dataframe, key = var3, value = value, var1, var2)

>

Task: Use gather() to create a new data frame fig1.4.data from share.data in which the different shares shall be the source columns for a new key-value pair. Specify share_type as key. The associated values shall be stored in the column shares.

#< task
# Load the required package 'tidyr' from the library:
library(tidyr)
#>
fig1.4.data = gather(share.data, key = share_type, value = shares, -exempt_group)
#< hint
display("Have a look at the info-box on the 'gather()' function.")
#>

Task: This data set can now be used to generate a bar graph with the different exemption categories exempt_group on the x-axis and the shares on the y-axis. Fill in the missing parts of the code and generate the graph.

#< task_notest
# Setting 'fill = share_type' colors the different share types under 
# 'share_type' differently. By specifying 'stat = identity' the hight 
# of the bar represents the values in the data (instead of counting 
# the number of cases in each group (default)). 'position = dodge' seperates 
# the different 'fill-categories' in different bars (instead of stacking 
# them on top of each other (default)).

# fig1.4 = ggplot(data = ??? , aes(x = ??? , y = ??? , fill = share_type)) +
#         geom_bar(stat = "identity", position = "dodge", color = "Black") +
#         scale_fill_manual(values=cbPalette)
# 
# fig1.4
#>
fig1.4 = ggplot(fig1.4.data, aes(x=exempt_group, y=shares, fill=share_type))+
        geom_bar(stat = "identity", position = "dodge", color = "Black") +
        scale_fill_manual(values=cbPalette)

fig1.4

Figure 1.4: Relative shares in exemption categories. Source: Adapted from Martin et al. (2014, p. 2489).

Figure 1.4 illustrates the shares the five different exemption categories take in terms of number of firms (blue), employees (yellow) and $CO_2$ emissions (grey).

Give the share of $CO_2$ emissions of EU ETS firms (in %) that is actually subject to full auctioning after the Carbon Leakage Decision:

< quiz "Shares in exemption categories"

question: answer: 16 roundto: 1

>

In the description of Figure 1.4 we will concentrate on the most important aspects. We can see that the share of $CO_2$ emissions subject to full auctioning amounts to just 16%. The remaining emissions are exempt by the EC's criteria. Category $B\: and\: CI>5\;\%$ comprises with 25% the largest share of exempted $CO_2$ emissions. The largest amount of jobs affected by the exemptions belongs to sectors in category $B\: and \: CI<5\;\%$. Category $B$ as a whole takes the largest share in terms of number of firms, employees and $CO_2$ emissions compared to the other exempted groups. Conversely, category $A$ takes the smallest share in every aspect compared to categories B and C. In consequence more firms and emissions are exempt from full auctioning due to trade intensity (category $B$) than to carbon intensity (category $A$), and these firms comprise the largest amount of jobs as well.

< info "Related Literature"

The share of $CO_2$ emissions that is not exempt from full auctioning is not clear cut, but depends on the measurement and calculation method. The 16% of emissions that are subject to full auctioning we derived above can be seen as the lower bound of figures found in the related literature. Juergens, Barreiro-Hurlé and Vasa (2013) found a share of 23% in their assessment. They used estimations of direct emission figures, as not all emission data is accessible due to confidentiality restrictions in CITL.

>

iii) Summary and discussion

There are two key insights we should take away from this exercise. First, by giving in to the threat of carbon leakage the EC has established very generous exemption criteria, exposing only 16% of $CO_2$ emissions to full auctioning. In consequence the ownership of emissions is left with European industries and is not transferred to taxpayers. The polluter-pays principle underlying European environmental policy-making is not implemented in this case. Second, most of the sectors are exempt from auctioning on the basis of their high trade intensity. Martin et al. argued that the link between a credible relocation threat and the EC's definition of trade intensity is not clear-cut. The TI measure is defined as the sum of import and export penetration (see info-box). The use of import penetration as a proxy for competitive disadvantages is justifiable: If imports from non-European countries increase, EU products are consequently substituted by relatively cheaper non-EU products. The additional cost imposed by the carbon policy cannot be passed through to the customers. Firms might relocate due to this competitive disadvantage. But the TI figure contains export penetration as well. This measure strongly depends on country specific factors like natural resource deposits or skilled labour force. These factor specificities create an absolute advantage that will not be diminished by higher carbon cost. The producer of Swiss watches for instance cannot relocate to China. A high TI measure may be due to a high export penetration rather than a high import penetration. In this case it might not be necessary to exclude firms from permit auctioning since no credible relocation threat arises.

These insights are our motivation to further investigate the free-permit allocation schemes established under the Carbon Leakage Decision and compare them to Martin et al.'s alternative proposal.

To move to Exercise 2 you can click on the Go to next exercise...-button.

The graphs and the associated discussions developed in this exercise can be found on pages 2486-2490 of the paper.

! start_note "Data Preparation Exercise 1"

The code chunk below shows the code to create the data frame firm.data used in Exercise 1 from the original data frame basicdata.dta which was provided by Martin et al.

#< task_notest
# To read files in Stata version 12 into a data frame we need the 
# 'foreign' package. 
library(foreign)
basic.data = read.dta("data/basicdata.dta")

# Create 'firm.data' by extracting only manufacturing firms that are
# part of the EU ETS and for which sectoral information are available.
# Reduce the data frame to variables that are used in Exercise 1 and
# assign more intuitive variable names where necessary.
firm.data = basic.data %>%
  filter(!(is.na(sec4dig)),
         is.na(nonmanufacturing),
         is.na(notETS)) %>%
  select(sec4dig,
         carb_int = vv_xxx0,
         trade_int = tt_xxx0,
         n_install = ninstallations,
         emp = empBigorbis,
         co2 = surr,
         firm_sec = countid)%>%
  mutate(trade_int = ifelse(is.na(trade_int),0,trade_int))
#>

! end_note

Exercise 2 The Vulnerability Score

As we have seen in Exercise 1, Martin et al. questioned that the exemption criteria proposed under the EU ETS are related to actual relocation risk. They therefore developed their own vulnerability score to measure a firm's vulnerability to carbon pricing. The improved compensation rules we are going to develop in Exercise 3 will be based on this score. This exercise will introduce you to the key dimensions of this vulnerability score. In the end of the exercise the score will be used to test the validity of the EC's exemption criteria.

i) Translation of interview answers into a score

In order to ascertain how European industries react to carbon pricing and climate policy measures in general, Martin et al. conducted telephone interviews with managers of 761 manufacturing firms. Interviews were carried out in six European countries: Belgium, France, Germany, Hungary, Poland and the United Kingdom. To obtain information on the critical subject of a firms downsizing decision in response to carbon pricing, interviewers made use of a survey tool developed by Bloom and van Reenen (2010). One feature of this survey tool is to ask open-end-questions and let the interviewers assign scores to the answers given. In combination with other elements of Bloom and van Reenen's survey format, this helps to avoid common bias sources in interviews. To determine a firm's expected relocation or outsourcing behaviour due to carbon pricing, researchers asked:

Do you expect that government efforts to put a price on carbon emissions will force you to outsource part of the production of this business site in the forseable future, or to close down completely? (Martin et al., 2014, p. 2491)

Interviewers translated the answers to this question into the socalled vulnerability score (VS). The score is defined on a scale from 1 to 5 and is assigned as follows:

Asignment of interview answers to the VS. Figure 2.1: VS mapping. Source: Author's own graph.

In the following, we want to analyze how managers responded to this question. Therefore, we need to load data containing the interview results.

Task: Use the read.table() command to load the data frame vs.data and assign it to a new variable vs.data. Remember that data frames are stored in the directory in a folder called data. Display the first entries of vs.data using the head() command.

vs.data = read.table("data/vs.data")
head(vs.data)
#< hint
display("The basic construction of the commands is 'vs.data = read.table('name of data frame')' and 'head(vs.data)'.")
#>

Each row contains the information on one firm. Besides the VS, the data set contains further firm-level information that was matched to the VS from other databases we already encountered in Exercise 1.

< info "vs.data"

The data frame contains information on 761 manufacturing firms in six European countries.

The data preparation code can be found in the footnote of this exercise.

>

ii) Mapping of the VS

A closer look on the chosen mapping from managers' relocation expectations into the VS, illustrated in Figure 2.1, reveals that small relocation probabilities are covered by the score in more detail than high relocation probabilities. Namely, the first two fifths of the score cover relocation probabilities from zero to 10%, while the other three fifths of the score cover the remaining 90%. This concentration on low relocation probabilities is justifiable if overall relocation probabilities in response to carbon pricing are low. Let us therefore have a look at the distribution of the VS in a histogram.

Task: Hit the check-button to generate a graph that shows the relative frequency of different score realizations.

#< task
# 'geom_histogram()' creates a histogram, 'binwidth = 1' ensures that 
# each score realization is collected in one bin,
# '(..count..)/sum(..count..)' displays relative frequencies
fig2.2 = ggplot(data = vs.data, aes(x=vs)) +
  geom_histogram(binwidth = 1, aes(y = (..count..)/sum(..count..))) +
  ylab("Share of firms") + 
  xlab("Vulnerability Score")
fig2.2
#>

Figure 2.2: VS distribution. Source: Author's own graph.

The histogram in Figure 2.2 shows that most managers expect pricing of carbon emissions to have very little impact on the relocation behaviour of their company. Over 60% believe that carbon pricing will have no detrimental effect at all. A score of three or higher was assigned to the answers of less than 29% of managers, meaning that less than a third of firms is expected to relocate more than 10% of their business.

The VS attempts to measure specifically the relocation tendency of firms in response to carbon pricing. But relocation decisions depend on many more aspects than only $CO_2$ emission costs. Examples of other important factors might be the stability of investment conditions or the availability of skilled labour (European Commission, 2016). Therefore, it is plausible that we observe a high relative frequency of low VS.

In later exercises we will examine how the VS of a firm changes when it can expect to receive free pollution allowances for 80% of its $CO_2$ emissions. To be able to measure a gradient of the score, it is important that the score can differentiate changes with sufficient precision.

On these grounds it is reasonable that Martin et al. tried to capture low relocation probabilities in more detail than high ones.

iii) Descriptive and test statistics of the VS

Let us have a look at further descriptive statistics of this VS.

Task: Calculate the sample mean and standard deviation of the VS across all firms. You can use the empty code chunk below to perform your calculations. Take into account that for some firms in the sample the VS might not be available. Pressing the hint-button provides a solution to this task. State your answers in the quiz below!

#< task_notest

#>
#< hint
display("Execute 'mean(vs.data$vs, na.rm = T)' and 'sd(vs.data$vs, na.rm = T)'")
#>

< quiz "Sample mean and standard deviation of the VS across all firms"

parts: - question: 1. State the sample mean up to 2 digits answer: 1.87 roundto: 0.01 - question: 2. State the sample standard deviation up to 2 digits answer: 1.29 roundto: 0.01

>

As already expected after assessing the histogram, the mean VS across all firms is rather low, namely 1.87 with a standard deviation of 1.29.

You may have noted that the data frame vs.data contains a variable ets. This variable indicates whether the interviewed firm takes part in the EU ETS or not. Martin et al. conducted interviews with both groups. In a next step, we want to know if firm managers have different expectations on the relocation behavior of their company depending on their EU ETS participation status.

Task: Use the commands summarize() and group_by to calculate the VS group means of ETS and non-ETS firms. 

vs.data %>% 
  group_by(ets) %>% 
  summarize(mean(vs, na.rm = TRUE))
#< hint
display("As introduced in Exercise 1, a possible solution might have the following form: 'datafram %>% group_by(var1) %>% summarize(mean(var2, na.rm = TRUE))'")
#>

Notably, the mean VS of non-ETS firms is lower than that of ETS firms (1.49 versus 2.15). In a next step, we want to know if this observed difference in means is statistically significant. If certain conditions are met, a two sample t-test is the most powerful tool to test the null hypothesis stating that the sample means of the two groups are equal. One of these conditions requires that the two samples being compared follow a normal distribution. Figure 2.2 reveals that the normality assumption is most likely violated by the samples. To avoid misleading results, we will therefore conduct a non-parametric Wilcoxon-Mann-Whitney test. This test examines whether the sample distributions of the two groups are identical (Sheskin, 2007, chapter Test 12). Thereby, it does not make assumptions about the specific underlying distribution, especially it does not have to be normal. It tests the following null hypothesis:

$H_0$: The distributions of VS for the two groups are identical.

Task: The code to conduct a Wilcoxon-Mann-Whitney test is already provided. Just click check to execute the test.

#< task
# Wilcoxon-Mann-Whitney test: 
wilcox.test(formula = vs~ets, data = vs.data)
#>

The Wilcoxon-Mann-Whitney test reports a p-value of $p=1.5\cdot 10^{-9}$. Thus, we can confidently dismiss the null hypothesis that the distributions of VS of ETS and non-ETS firms are identical and conclude that there is true difference in the mean VS. As one might expect, ETS firm managers assume more detrimental effects caused by carbon pricing than those of the reference group. In the optimizations we are going to conduct in Exercises 3 and 4, we are therefore going to concentrate on ETS firms only.

Since the EU ETS is a European trading system, we are interested in managers' expectations concerning carbon pricing across different countries. A very nice tool to receive meaningful information on different groups in a sample is the pirateplot() function in the yarrr package.

< info "pirateplot()"

The pirateplot() function enables us to obtain detailed but easy to grasp information by writing very little code. Its output is a so called RDI plot (Raw data, Description and Inference). These plots combine features of both box plots and bean plots. For an example, have a look at the next task. Further information are provided in the book YaRrr! The Pirate's Guide to R.

>

Task: Have a look at the provided example. It creates a plot showing the distribution of the VS depending on the country where it was obtained. Hit the check-button to create a first pirateplot. 

#< task
# Load the required package 'yarrr' from the library
library("yarrr")

pirateplot(formula = vs ~ country,
           data = vs.data,
           inf = "ci",
           xlab = "Country",
           ylab = "Vulnerability Score",
           main = "Distribution of VS across countries") +
  abline(h=mean(vs.data$vs, na.rm = T))
#>

Figure 2.3: Vulnerability to carbon pricing in different countries. Source: Adapted from Martin et al. (2014, p. 2492).

Note: Martin et al. used a regression analysis controlling for interview noise to determine the confidence intervals in their graph. They yielded slightly different results.

In Figure 2.3, the horizontal lines indicate the mean with the surrounding boxes illustrating the 95% confidence intervals. In the background we observe the raw data points. The beans represent the smoothed densities of the five raw score values obtained in the survey. As a reference the overall mean VS is included as a black horizontal line in the graph.

< quiz "Vulnerability score across countries"

question: Check the three countries which exhibit the highest VS. mc: - Belgium - France - Germany - Hungary - Poland* - United Kingdom success: Great, all answers are correct! failure: Not all answers correct. Try again.

>

It is apparent that the mean VS of all countries is lower than 3. Moreover, none of the 95% confidence intervals includes the score 3. This allows us to conclude that the risk of relocation due to carbon pricing is limited to less than 10% production outsourcing in all of the countries. The beans unveil that the distribution of scores is heavily skewed towards a score of 1 in all countries.

In a next step, we want to do a similar analysis concerning managers' relocation expectations in different industries.

Task: Complete the code in the chunk below to create a second pirate plot, which displays the distribution of VS depending on the industry (ets_sector) where it was obtained.

#< task_notest
# This command ensures that the entire graph is displayed. 
# You don't have to worry about it.
par(las = 2, mar = c(9.3,4.2,3,0))

# Enter your code below
# pirateplot(formula = ???,
#            data = vs.data,
#            inf = "ci",
#            xlab = "",
#            ylab = "Vulnerability Score",
#            main = "Distribution of VS across sectors") +
#   abline(h=mean(vs.data$vs, na.rm = T))
#>

pirateplot(formula = vs ~ ets_sector,
           data = vs.data,
           inf = "ci",
           xlab = "",
           ylab = "Vulnerability Score",
           main = "Distribution of VS across sectors") +
  abline(h=mean(vs.data$vs, na.rm = T))

Figure 2.4: Vulnerability to carbon pricing in different industries. Source: Adapted from Martin et al. (2014, p.2492).

Note: Martin et al. used a regression analysis controlling for interview noise to determine the confidence intervals in their graph. They yielded slightly different results.

< quiz "Vulnerability score across industries"

question: Check all sectors that have an average score above 3. mc: - Cement - Ceramics - Fuels - Glass - Other Minerals* - Vehicles success: Great, all answers are correct! failure: Not all answers correct. Try again.

>

Figure 2.4 reveals that none of the sectors displays a severe vulnerability to carbon pricing. Fuels, Glass, Iron and Steel, and Other Minerals are the industries with the highest VS. Thereby, Other Minerals displays the highest score with a mean VS greater than 3 (3.4). Only in the cases of Ceramics, Fuels, Glass, and Iron and Steel the value 3 is included in the 95% confidence interval. We can therefore not rule out the risk that firms in these industries shift at least 10% of their production outside the EU. In all other industries the average VS is substantially lower than 3. In no industry can we observe that the 95% confidence interval of average scores include complete relocation. We may therefore conclude that chances are low that an average firm closes down completely and moves to an unregulated country as a consequence of carbon pricing.

Overall, the industry a firm belongs to has a higher influence on the vulnerability to carbon pricing than the the country of origin. The country means range between 1.5 (Hungary) and 2.1 (Germany) and the differences between the countries appear to be rather minor, whereas industry means range between 1.26 (Machinery & Optics) (neglecting Construction since there are only three firms in the sample) and 3.4 (Other Minerals). Therefore, there is a higher variation in vulnerability between sectors than between countries.

The raw data and the smoothed densities reveal that in the identified most vulnerable industries the VS are rather evenly distributed. Consequently firms in Fuels, Glass, Iron and Steel and Other Minerals can be considered as heterogeneous with respect to their vulnerability to carbon pricing. Therefore, the EU's approach of basing exemption rules on sector specific measures like carbon or trade intensity may be inefficient.

iii) How does the VS relate to carbon and trade intensity?

The VS introduced above claims to be a direct measure of a firm's vulnerability to carbon pricing. Martin et al. (2014b) used it therefore in an accompanying paper as a tool to test the validity and accuracy of the EC's exemption criteria established in the Carbon Leakage Decision. Their basic idea was that TI and CI should be positively correlated with the VS to be good proxies for a firm's relocation risk. They used a variant of the following regression equation to test this hypothesis:

$$ VS_{i,s,c} = \beta_0 + \beta_{TI}\cdot TI_s + \beta_{CI}\cdot CI_s + \delta_c + \varepsilon_{i,s,c}.\qquad (2.1) $$

In this formula $VS_{i,s,c}$ is firm $i$'s vulnerability score, $TI_s$ and $CI_s$ are the sectoral trade and carbon intensity as introduced in Exercise 1. The vector $\delta_c$ contains a full set of country dummies to control for country effects. In addition to this basic set up Martin et al. (2014b) control for interview noise. Since the code for this regression is not publicly available and we therefore do not know which interview specific variables are included, we, unlike Martin et al., neglect these effects in our analysis.

It is important to emphasize that we do not intend to conduct a classical econometric analysis as focused on in many econometric textbooks. Especially we do not assume that equation (2.1) represents a model for the true underlying data generation process. Our aim is not to interpret the relationship between the VS and the CI and TI measures in a causal sense. We simply want to add another piece to the descriptive analysis of the VS, and gain information on its partial correlations with CI and TI while controlling for country effects.

< info "lm()"

The lm() function uses the method of ordinary least squares (OLS) to fit a linear relationship to data. OLS estimates the coefficients by minimizing the sum of squared residuals. The residuals are the differences between the actual data points and the estimated curve.

Let us assume that we have a data frame dataframe that contains the three different variables y, x1 and x2. To fit a linear model of the relationship between variable y and variables x1 and x2 we can apply the lm() function as follows:

lm(y ~ x1 + x2, dataframe)

>

Task: Click the check-button to use the lm() function to conduct the linear approximation of the VS of ETS firms displayed in equation (2.1).

#< task
reg = lm(formula = vs ~ carb_int + trade_int + country, 
          data = filter(vs.data, ets == 1))
#>

Before we visualize the results of this regression, let us reconsider our observations in the pirateplot displaying the distribution of the VS across sectors. The smoothed densities of the different sectors are differently shaped. Some, as Machinery & Optics, are heavily skewed towards a score of 1, others, as Other Minerals are rather evenly distributed. This is a strong indication for heteroskedasticity. Martin et al. (2014b) therefore used robust standard errors clustered by four digit sectors.

We are going to use the showreg() function from the regtools package to display the results of our regression. This function creates an elegant presentation of regression results and allows to calculate and display robust standard errors in an easy fashion.

Task: Click the check-button to calculate robust standard errors clustered by four digit sectors and to display the regression results using the showreg() function.

#< task
# Load the required package 'regtools' from the library
library(regtools)

# To calculate clustered standard errors 'showreg()' requires 
# the original data without missing values 
reg$custom.data = filter(vs.data, ets==1, !is.na(vs), !is.na(trade_int))

# Display regression results with robust standard errors 
# clustered by four digit sectors
showreg(list(reg), robust = T, robust.type = "cluster", 
        cluster1 = "sec4dig", digits = 3)
#>

The table above shows the estimated coefficients for all the variables of the right hand side of equation (2.1). Underneath the coefficient estimates we find the robust standard errors in parentheses. The asterisks behind the estimates indicate their statistical significance. Focussing on the estimated coefficients of TI and CI, we observe that in this linear approximation only CI is significantly positively associated with the VS. However, there is no statistically significant association between vulnerability and TI. Controlling for interview noise Martin et al. (2014b) found qualitatively similar results. The coefficient of determination $R^2 = 0.08$ shows that only 8% of the variance in the VS can be explained by this linear approximation.

A practical difficulty in analysing regression results is that different explanatory variables are naturally measured in different units. In regression (2.1) for instance CI is the ratio of carbon cost to gross value added, whereas TI is the ratio of the sum of exports to and imports from third countries to the total market size. Effects on the dependent variable $VS_{i,s,c}$ induced by changes in the explanatory variables $TI_s$ and $CI_s$ are therefore difficult to compare. The function effectplot() of the regtools package aims at a normalization of these effects to make effect sizes comparable. This normalization is realized by determining the effect on the dependent variable when an explanatory variable changes ceteris paribus from its 10% quantile to its 90% quantile.

Task: Click check to display the regression results in an effectplot.

#< task
# The 'effectplot()' function cannot deal with factor variables in the 
# regression formula. The 'expanded.regression.data()' command therefore 
# creates a new data frame in which the factors of 'country' are 
# converted into dummy variables.
effect.data = expanded.regression.data(vs ~ carb_int + trade_int + country,
                                       data = filter(vs.data, ets == 1))

# Rerun regression (2.1) explicitly using dummy variables for the countries.
effect.reg = lm(regression.formula(effect.data), data=effect.data)

# Display the effectplot
effectplot(effect.reg)
#>

Figure 2.5: Effects of explanatory variables on the VS. Source: Author's own graph.

The effectplot in Figure 2.5 displays the effect on the VS as the explanatory variables change form their 10% quantile to their 90% quantile. In the case of the country dummies the graph shows the effect on the VS as the dummies change from zero to one. Differences in the sign of the effect can be distinct by the coloring of the bar.

< quiz "Effectplot"

question: How much larger is the VS of a firm with a high CI (90% quantile) compared to a firm with a low CI (10% quantile)? answer: 0.32 roundto: 0.01

>

Figure 2.5 reveals that changes in both CI and TI, are relatively weakly associated with changes in a firm's vulnerability to carbon pricing. A firm with a high CI (90% quantile) displays on average a VS that is increased by only 0.32 compared to a firm with a low CI (10% quantile). The TI shows with an increase of 0.18 in the VS an even weaker effect. With the results of the pirateplots above in mind, it is worth considering that the country of origin of a firm shows for some countries (Germany (0.61) and Poland (0.43)) stronger effects on the VS than CI and TI.

To encounter the argument that a linear combination of continuous variables TI and CI approximating the VS is not suitable to evaluate the fixed exemption thresholds established by the EC, Martin et al. (2014b) adopted a second approach: They conducted a regression with the VS as the dependent variable and dummy variables representing the EC's exemption categories ($A$, $B$, $C$, $not\: exempt$) as introduced in Exercise 1 as the explanatory variables. The result of this second approach yielded similar findings as found above: Only very carbon intensive firms (category $A$) exhibit a mean VS significantly higher than firms not exempt from full auctioning (reference group).

These findings support the implications concluded from economic reasoning in Exercise 1: Especially trade intensity is not a good indicator for relocation risk.

iv) Summary and discussion

The findings above reveal that firms in all countries and most industries will relocate less than 10% of their business to less regulated countries in response to a price tag on carbon emissions. We may therefore conclude that the EC's policy of exempting 84% of emissions from full auctioning is not justifiable.

The graphical analysis of the different industries demonstrates that there is a great heterogeneity across sectors and across firms in vulnerability to carbon pricing. An approach that yields at exempting entire sectors as chosen by the EC might therefore be inappropriate.

Our regression analysis shows that trade intensity, in contrast to carbon intensity, is a bad proxy for a firms vulnerability to carbon pricing. This problem is especially severe since most firms (27%) are exempt from full permit auctioning on the basis of the TI measure.

The mapping and the descriptive statistics developed in this exercise can be found on pages 2491-2492 of the paper.

! start_note "Data Preparation Exercise 2"

The code chunk below shows the code to create the data frame vs.data used in Exercise 2 from the original data frame basicdata.dta which was provided by Martin et al.

#< task_notest
# Load the data frame 'basicdata.dta'
basic.data = read.dta("data/basicdata.dta")

# Create 'vs.data' by extracting one interview result per firm 
# (some firms were interviewed more than once).
# Reduce data frame to variables that are used in Exercise 2 and
# assign more intuitive variable names where necessary.
vs.data = basic.data %>% 
  filter(unique ==1)%>%
  select(vs = fimpact_score_clean,
         country,
         ets_sector = mcetsdig,
         sec4dig,
         carb_int = vv_xxx0,
         trade_int = tt_xxx0,
         ets = ETS_id)
#>

! end_note

Exercise 3.1 - Optimal Allocation - A Model

In Exercise 1 we discussed that by putting a price tag on $CO_2$ emissions, the EU faced the dilemma of introducing an effective environmental policy on the one hand and protecting domestic industries against the detrimental economic implications of this policy on the other hand. If firms decided to relocate to non-regulated countries in response to these detrimental implications, this would cause damage to both the overall goal to reduce carbon emissions and the domestic economy. In the course of this problem set we will follow Martin et al.'s approach and focus on two concrete variants of relocation damage: Where applicable, we will differentiate between the amount of jobs lost in the EU and the amount of leaked $CO_2$ emissions caused by relocation. As seen in Exercise 1, the EC decided to allocate a fixed amount of pollution permits to firms for free to encounter the relocation threat. In this exercise we want to develop Martin et al.'s theoretical framework on how to distribute these free permits in order to minimize the aggregate expected damage of relocation.

i) A single firm's contribution to aggregate relocation risk

Martin et al. based their model on the relocation decisions of individual firms located in an ETS-country. They assumed that these relocation decisions depend on the profit impact of carbon pricing. Therefore, they described the profit $\pi_i (p,q_i)$ of such a firm $i$ as a function of the permit price $p$ and the number of free permits $q_i$ it receives. Since the EU ETS is a cap-and-trade system we can consider the overall amount of pollution permits as exogenously fixed. The first person to develop a theoretical model for permit prices under these premises was Montgomery (1972). He found that permit trading reduces emissions cost-optimally and that permit prices are equal to firms' marginal abatement costs in the optimum. In the course of their model Martin et al. assumed the overall distribution of firms' abatement costs to be constant. As a consequence the permit price $p$ does not change and can be ignored in the profit function for better legibility: $\pi_i (q_i)$.

Do you expect the partial derivative of a firm's profit with respect to the number of free permits $\frac{\partial \pi_i (q_i)}{\partial q_i}$ to be positive or negative?

< quiz "Profit"

question: sc: - positive* - negative success: Great, your answer is correct! failure: Try again.

>

Martin et al. assumed a positive additional profit if a firm receives additional free permits. Since permits are a scarce resource and need to be acquired to emit $CO_2$, free permits can be considered as a lump-sum subsidy to firms.

In this model a firm $i$ will relocate if its profit $\pi_i$ in an ETS-country is smaller than the difference of profits $\pi_{if}$ it could earn in an unregulated country $f$ and incurring relocation costs $\kappa_i$: $\pi_i (q_i) < \pi_{if} - \kappa_i$. Martin et al. assumed that while regulators can observe a firm's profit in its home country, they have no detailed information on its net relocation costs $\varepsilon_i \equiv \kappa_i - \pi_{if}$. In their model regulators only have the knowledge that $\varepsilon_i$ is an independent and identically distributed random variable. They can observe its mean $\mu_{\varepsilon}$ and standard deviation $\sigma_{\varepsilon}$ and know that it follows the differentiable cumulative distribution function $F_i(\cdot)$. Subsequent to these assumptions we can define the binary relocation variable $$ y_i \equiv \mathbb{1}{\varepsilon_i < -\pi_i(q_i)}, $$ that takes the value 1 if the condition is fulfilled and firm $i$ relocates and 0 otherwise. Given firm $i$ receives $q_i$ permits for free, its relocation probability can be defined as $\mathrm{Pr}(y_i =1 \mid q_i) = F_i[-\pi_i(q_i)]$. Since we are interested in minimizing the aggregate expected damage of relocation, the relocation probability is weighted by the damage $damage_i$ which the relocation of firm $i$ causes. Damage is either expressed in terms of number of employees $emp_i$ or tons of $CO_2$ emissions $co2_i$ at firm $i$. Firm $i$'s contribution to aggregate relocation risk is therefore given by $$ risk_i (q_i) = F_i[-\pi_i(q_i)]\cdot damage_i, \quad damage_i \in [emp_i, co2_i].\qquad (3.1.1) $$ Depending on which objective we want to minimize - Carbon Leakage Risk or Job Risk - $damage_i$ takes either the values $co2_i$ or $emp_i$. The above formulation implies that if a firm decides to relocate, "all of its jobs are lost and all of its emissions 'leak' to nonregulated countries" (Martin et al., 2014, p.2498).

ii) Minimizing aggregate relocation risk

In this framework the EC's optimization problem is to minimize the sum of the individual firm's relocation risks. The overall amount of free permits is constrained by the EC to $\bar{Q}$. Formally we can write: $$ \min_{q_i \geq 0} \sum_{i=1}^n risk_i (q_i) \qquad s.t. \sum_{i=1}^n q_i \leq \bar{Q}. \qquad (3.1.2) $$

Remember that $F_i[-\pi_i(q_i)]$ was assumed to be a continuously differentiable cumulative distribution function. Do you think its first derivative with respect to the number of free permits $\frac{\partial F_i[- \pi_i (q_i)]}{\partial q_i}$ is positive or negative?

< quiz "Marginal relocation probability"

question: sc: - positive - negative* success: Great, your answer is correct! failure: Try again.

>

Applying the chain rule reveals that an additional free permit always results in a marginal reduction in firm $i$'s relocation probability. In consequence, the permit constraint of (3.1.2) holds with equality. If you are interested in the proof have a look at the next info-box.

< info "Proof of binding constraint"

In the set up of the model we assume: $F_i$ is a continuously differentiable cumulative distribution function. Cumulative distribution functions have the general property to be monotonically increasing: $F_i(x_1)\leq F_i(x_2)\:\forall\: x_1 < x_2$. $\pi(q_i)$ is monotonically increasing. Thus, $\frac{\partial F_i[- \pi_i (q_i)]}{\partial q_i}= -F_i^{\prime} [-\pi_i(q_i)]\frac{\partial \pi_i(q_i)}{\partial q_i} < 0$, and $F_i$ is monotonically decreasing. * Damage caused by relocation is non-negative: $damage_i \geq 0$.

Therefore we can proof by contradiction that the constraint of (3.1.2) is binding in the optimum.

Proof:

Assume the constraint is not binding in the optimum. Let $q_1^,\cdots ,q_n^$ be a solution of the minimization and $$ \sum_{i=1}^n q_i^ < \bar{Q}. $$ Then, there exists a small variation $\Delta q >0$ such that $$ \sum_{i=1}^n (q_i^ +\Delta q) \leq \bar{Q} $$ and since $F_i$ is monotonically decreasing $$ \quad \sum_{i}^n F_i[-\pi_i(q_i^ +\Delta q)]\cdot damage_i < \sum_{i}^n F_i[-\pi_i(q_i^ )]\cdot damage_i. $$ It follows that $q_1^,\cdots ,q_n^$ cannot be a solution. qed

>

The shadow price $\lambda^*$ in this minimization problem tells us by how much the optimal value of aggregate relocation risk increases if the total number of free permits is decreased by one. Following our previous findings on the first derivative of the relocation probability function $F_i[-\pi_i(q_i)]$, we can state that the shadow price is positive. We can further write the first order condition for an interior solution as

$$ \lambda = F_i^{\prime} [-\pi_i(q_i)]\frac{\partial \pi_i(q_i)}{\partial q_i}\cdot damage_i \quad \forall i.\qquad (3.1.3) $$ For more details on the derivation of this expression (3.1.3) and shadow prices have a look at the info-box on shadow prices.

< info "Shadow prices"

Since we have shown that the constraint holds with equality in the optimum, we can use the method of Lagrangian multipliers to rewrite the optimization problem. The Lagrangian of this optimization can be defined as $$ \mathcal{L}(q_i , \lambda) = \sum_{i=1}^n risk_i (q_i) - \lambda\left(\bar{Q} - \sum_{i=1}^n q_i \right). $$ The first order conditions yield: $$ \frac{\partial\mathcal{L}(q_i , \lambda)}{\partial q_i} = risk_i^{\prime}(q_i) + \lambda = - F_i^{\prime} [-\pi_i(q_i)]\frac{\partial \pi_i(q_i)}{\partial q_i}\cdot damage_i +\lambda \overset{!}{=} 0 $$ $$ \frac{\partial\mathcal{L}(q_i , \lambda)}{\partial \lambda} = \sum_{i=1}^n q_i - \bar{Q} \overset{!}{=} 0 $$ Therefore, we can write $$ \lambda = risk_i^{\prime}(q_i) = F_i^{\prime} [-\pi_i(q_i)]\frac{\partial \pi_i(q_i)}{\partial q_i}\cdot damage_i. $$ In general, $\lambda$ measures the sensitivity of the Lagrangian to variations in the constraint parameter $\bar{Q}$: $$ \frac{\partial\mathcal{L}(q_i , \lambda)}{\partial \bar{Q}} = \lambda. $$ In the optimum, the Lagrangian multiplier allows an interesting interpretation. Assume we have a solution to our optimization problem consisting of the optimal value of the objective function $\sum_{i=1}^n risk_i (q_i^)$, the optimal variable values $q_1^,\cdots,q_n^$ and the resulting value of the Lagrangian multiplier $\lambda^$. Then it can be shown that the partial derivative of $\sum_{i=1}^n risk_i (q_i^)$ with respect to $\bar{Q}$ is equal to the resulting value of the Lagrangian multiplier $\lambda^$: $$ \frac{\partial\sum_{i=1}^n risk_i (q_i^)}{\partial \bar{Q}} = \lambda^. \qquad (3.1.3.1) $$ In our optimization $\lambda^$ therefore tells us by how much aggregate relocation risk changes if the total number of free permits is decreased by one. In economics, when the objective function displays a monetary value, $\lambda^$ is often referred to as the shadow price. In the case of our optimization, the price has no direct monetary interpretation but represents increased relocation risk.

See Intriligator (2002, chapter 3) for a more formal definition and description of the equations used above and a proof of equation (3.1.3.1).

>

The first order condition (3.1.3) implies that the regulator has to distribute free permits among firms such that the allocation equalizes marginal expected relocation damages between firms. In consequence free permits should be allocated to firms, where they generate the largest reduction in contribution to aggregate relocation damage. Martin et al. claimed to be the first to bring this simple approach up in the context of free permit allocation.

iii) Using the VS to specify the relocation probability function $F_i$

So far we have made very broad assumptions on the relocation probability function $F_i$ to develop the model. In order to apply the model to data, we need to become more specific. Martin et al. assumed firm $i$'s profit increases linearly in the amount of free permits it receives: $$ \pi_i (q_i) = \delta_{0i}+\delta_{1i}q_i.\qquad (3.1.4) $$ Further, they assume that the net cost of relocation $\varepsilon_i$ follows a logistic distribution. We can therefore specify firm $i$'s relocation probability to be

$$ \mathrm{Pr}(y_i = 1|q_i) = F_i (-\pi_i (q_i)) = \frac{1}{1+\exp(\beta_{0i}+\beta_{1i}q_i)}.\qquad (3.1.5) $$ For details see the info-box on the Logistic Distribution.

< info "Logistic Distribution"

The cumulative distribution function of the logistic distribution is given by $$ F(x;\mu,s) = \frac{1}{1+\exp \left(-\frac{x-\mu}{s} \right)},\quad -\infty < x < \infty, $$ with mean $\mu$ and scale parameter $s=\frac{\sqrt{3}\sigma}{\pi}$, which is proportional to the standard deviation $\sigma$ (Balakrishnan, 1992, page 2).

With equation (3.1.4) we receive $$ F_i (-\pi_i (q_i)) = \frac{1}{1+\exp \left(-\frac{-\delta_{0i}-\delta_{1i}q_i-\mu_{\varepsilon}}{s_{\varepsilon}} \right)} = \frac{1}{1+\exp \left(\frac{\delta_{0i}+\mu_{\varepsilon}}{s_{\varepsilon}}+\frac{\delta_{1i}}{s_{\varepsilon}}q_i \right)}, $$ and can define parameters $\beta_{0i}\equiv \frac{\delta_{i0}+\mu_{\varepsilon}}{s_{\varepsilon}}$ and $\beta_{1i} \equiv \frac{\delta_{1i}}{s_{\varepsilon}}$. As stated above, regulators have knowledge about the mean $\mu_{\varepsilon}$ and the standard deviation $\sigma_{\varepsilon}$ of firm $i$'s net cost of relocation $\varepsilon_i$.

>

In equation (3.1.5) above, we have two unknowns that need to be classified using the interview data: $\beta_{0i}$ and $\beta_{1i}$. Since there are two unknowns, we need two equations. In Exercise 2, the Vulnerability Score was introduced. We have seen how this score captures managers' expectations on relocation probabilities under the assumption that no free permits exist: $q_i =0$. We can use these relocation expectations to determine $\beta_{0i}$ as follows: Given firm $i$ receives zero free permits, its relocation probability becomes $$ \mathrm{Pr}(y_i = 1|q_i = 0) = \frac{1}{1+\exp(\beta_{0i})}. $$ Rearranging this equation yields $$ \longrightarrow\quad \beta_{0i} = \log \left(\frac{1-\mathrm{Pr}(y_i = 1|q_i = 0)}{\mathrm{Pr}(y_i = 1|q_i = 0)} \right). \qquad (3.1.6) $$ Have a look at the exact mapping between the VS and relocation probabilities as introduced in Exercise 2 in the info-box below.

< info "Mapping from relocation probabilities into the VS"

Managers' assessment of their companie's relocation probabilities are mapped as follows into the vulnerability score:

Mapping of managers' assessment of their companie's relocation probability into the VS. Thereby, relocation probabilities belonging to scores 1, 3 and 5 are assigned by definition. Relocation probabilities of intermediate scores 2 and 4 are interpolated.

>

Managers were asked as well how these expectations would change if firms were granted to receive free pollution permits for 80% of their emissions. Using the same mapping, answers are assigned to a Vulnerability Score again. With this second score and the results for $\beta_{0i}$ we can determine $\beta_{1i}$:

$$ \mathrm{Pr}(y_i = 1|q_i = 0.8\cdot co2_i) = \frac{1}{1+\exp(\beta_{0i} + \beta_{1i}\cdot 0.8\cdot co2_i)} $$ $$ \longrightarrow\quad \beta_{1i} =\frac{1}{0.8\cdot co2_i}\left( \log \left(\frac{1-\mathrm{Pr}(y_i = 1|q_i = 0.8\cdot co2_i)}{\mathrm{Pr}(y_i = 1|q_i = 0.8\cdot co2_i)} \right)-\beta_{0i}\right).\qquad (3.1.7) $$

Therefore, firm $i$'s contribution to aggregate relocation risk is fully determined and we can rewrite equation (3.1.1). Besides substituting equation (3.1.5) into equation (3.1.1) we formally change the argument of the $risk_i$ function from permits $q_i$ to batches $b_i$. This accounts for the fact that not single emission permits $q_i$ are allocated to the firms, but rather batches $b_i$ comprising many permits (we will assume that one batch contains 15,000 permits, which allows the permit holder to emit 15,000 tons of $CO_2$). This substitution does not severely change the nature of the optimization problem, since batches $b_i$ are a linear transformation of permits $q_i$: $$ risk_i (b_i) = F_i(b_i) \cdot damage_i = \frac{1}{1+\exp (\beta_{0i} + \beta_{1i}\cdot batchsize\cdot b_i)} \cdot damage_i. \qquad (3.1.8) $$

iv) Example and function design

In the following, we want to apply the equations we have developed to a real firm to make the theory more tangible. For that purpose I randomly selected one firm from our data base. The firm has the following characteristics:

Task: Use the provided data and equations (3.1.6) and (3.1.7) to calculate $\beta_0$ and $\beta_1$ for that firm in the next code chunk. Do not forget to display your results. 

#< task_notest
# beta0 = ???
# beta1 = ???
#>
beta0 = log((1-0.99)/0.99)
beta1 = (log((1-0.01)/0.01)-beta0)/(0.8*2038472)  
beta0 
beta1
#< hint
display("The Logarithm is a predevined function in R. To calculate the Natural Logarithm of a 'number' simply type 'log(number)'. The mapping between VS and relocation probabilities can be found in the respective info-box above.")
#>

Now, that we have calculated the parameters $\beta_0$ and $\beta_1$ for our example firm, we can specify its relocation probability function as $$ F_{ex}(b_{ex}) = \frac{1}{1+\exp (-4.595 + 5.635\cdot 10^{-6}\cdot batchsize\cdot b_{ex})}.\qquad (3.1.9) $$

In a next step, we want to determine the expected number of jobs at risk of relocation if the example firm receives 20 batches of free permits. These 20 batches cover roughly 15% of that firm's $CO_2$ emissions.

Throughout this problem set we have already used a lot of different functions. Some of them were part of base R, like the log() function we used above, others were written by third parties and stored in special packages, like the ggplot() function from the ggplot2 package we used extensively in the last exercises. In R you have however the possibility to write your own functions as well. This is especially handy if non of the predefined functions suits your requirements and you want to do the same calculation or data manipulation multiple times. Since we will calculate the contribution to aggregate relocation risk for many different firms and batch allocations, it is a good idea to write a function which evaluates equation (3.1.8).

< info "Functions in R"

Functions in R generally have three components (Wickham, 2015):

The function returns the last value generated within the function. One can however use a return() statement anywhere in the function, which specifies the value to be returned.

>

Since all the information on firms we want to analyze using the model above are stored in data frames, our self-written function needs to adapt to the structure of these data frames. To apply this function to our example firm, its characteristics need to be stored in a similarly designed data frame.

Task: Hit the check-button to create a new data frame firm.ex using the data_frame() function. 

#< task
firm.ex = data_frame(beta0 = beta0,
                    beta1 = beta1,
                    emp = 667,
                    co2 = 2038472)
firm.ex
#>

The data frame firm.ex contains all the necessary information on our example firm to evaluate equation (3.1.9).

Task: Have a look at the next code chunk to see an example of how to define equation (3.1.8) as the function risk(). Hit the check-button to implement the function.

#< task
risk = function(data, b, damage){
  batchsize = 15000
  return(data[,damage]/(1+exp(data[,"beta0"]+data[,"beta1"]*batchsize*b)))
}
#>

The risk() function takes three arguments: A data frame data that contains all the relevant firm information, the number of batches b that is allocated to a firm, and a string damage that specifies the damage weight of interest ("emp", "co2"). The expected relocation risk is returned.

Task: Applying the risk() function, calculate the expected number of jobs at risk of relocation if the firm receives 20 batches of free permits. 

risk(data = firm.ex, b = 20, damage = "emp")
#< hint
display("Call the function 'risk(data = ???, b = ???, damage = ???)' and make the necessary assignments.")
#>

The result reveals that if this firm receives free permits covering 15% of its $CO_2$ emissions, it can be expected that 632 of its 667 total jobs are at risk to be relocated to an unregulated country.

v) Cost minimization

Conversely to the previous optimization aim, the EC could minimize the number of free permits allocated to firms by constraining the relocation risk to a certain level $\bar{R}$. Formally we can write: $$ \min_{q_i \geq 0} \sum_{i=1}^n q_i \qquad s.t. \sum_{i=1}^n risk_i(q_i) \leq \bar{R} \qquad (3.1.10) $$

Instead of minimizing total relocation risk by optimally allocating permits of a fixed amount of total free permits, program (3.1.10) minimizes total permits by optimally allocating risk of a fixed amount of total relocation risk. Assuming auctioning to be the basic principle of allocation, reducing the number of free permits decreases foregone auction revenue. Therefore, Martin et al. referred to this set up as the taxpayer's cost minimization.

Martin et al. found that, as for the primal program, the first order condition to optimization (3.1.10) yields that the marginal impact on relocation risk needs to be equalized across all firms receiving non-zero allocations by the last free permit.

Since equation (3.1.5) is strictly monotonic in $q_i$ we can invert equation (3.1.1) and get $$ q_i(risk_i) = \frac{1}{\beta_{1i}}\left( \log \left( \frac{damage_i}{risk_i}-1\right)-\beta_0\right) \quad if \quad risk_i < \frac{damage_i}{1+\exp (\beta_{0i})} \qquad (3.1.11) $$ $$ q_i(risk_i) = 0 \quad if \quad risk_i \geq \frac{damage_i}{1+\exp (\beta_{0i})}. $$ This equation can be understood as follows: Depending on the amount of $risk_i$ that is allocated to a firm $i$ this firm adds $q_i$ free permits to the total amount free permits.

Equivalently to aggregating permits into batches in the primal program, we aggregate single jobs or tons of $CO_2$ emissions at risk into predefined risk pieces of a certain size. Depending on the constraining risk figure, job risk or carbon leakage risk, risk pieces take different sizes. Martin et al. assumed that one piece of job risk contains 5 jobs, whereas one piece of carbon leakage risk contains 2000 tons of $CO_2$ emissions. In order to be used in a numerical solution procedure we need to write equation (3.1.11) as a function permits() similarly to the risk() function we have developed above. This is done in the next code chunk.

Task: Hit the check-button to implement equation (3.1.11) as the function permits().

#< task_notest
permits = function(data, pieces, damage){
 piecesize = ifelse(damage == "co2",2000,5)
 return(
 ifelse(piecesize*pieces<data[,damage]/(1+exp(data[,"beta0"])),
   1/data[,"beta1"]*(log(data[,damage]/(piecesize*pieces)-1)-data[,"beta0"]),
   0))
}

#>

The permits() function takes three arguments: A data frame data that contains all the relevant firm information, the number of risk pieces pieces that is allocated to a firm, and a string damage that specifies the damage weight of interest (“emp”, “co2”). The expected amount of permits necessary to achieve the risk level defined by pieces is returned.

To make this more tangible reconsider our example firm from above.

Task: Calculate the amount of free permits our example firm should receive in order to relocate 50 jobs (10 pieces of job risk). Use the permits() function defined above.

permits(data = firm.ex, pieces = 10, damage = "emp")
#< hint
display("Call the function 'permits(data = ???, pieces = ???, damage = ???)' and make the necessary assignments.")
#>

The result reveals that the firm under consideration needs to receive 1,261,285 free permits (not permit batches) to relocate 50 jobs to an unregulated country.

The model developed in this exercise can be found on pages 2497-2500 of the paper.

Exercise 3.2 Optimal Allocation - Dynamic Programming

In this exercise we want to develop an algorithm that solves our core problem of distributing free permits among firms with the objective to minimize aggregate relocation risk, as defined in Exercise 3.1. The attempt uses dynamic programming techniques. We start by looking at a small example consisting of only 5 batches of free permits and three different firms. The objective of this example will be the minimization of total job risk. Therefore, we need to figure out how many batches (if any) to allocate to each of these firms. Solving this example step by step will leave us with a better intuition on how the dynamic programming strategy works. Later on this procedure is put together to one generalized program which can be applied to all firms and possible permit allocations in our data set.

Note: If you are primarily interested in the economic content of the paper, you might want to skip this exercise which focuses on technical aspects of dynamic programming.

i) A single firm's contribution to aggregat relocation risk

As stated above, we will discuss a minimal example comprising three firms and 5 batches of free permits. Each batch consists of 15,000 permits, which is equivalent to the right to emit 15,000 tons of $CO_2$. Batches cannot be further subdivided, which means that the number of batches allocated to each firm is an integer.

Task: Assign the information given in this paragraph to the variables displayed in the window below.

#< task_notest
#no.firms = ???
#no.batches = ???
#> 
no.firms = 3
no.batches = 5

In Exercise 3.1 we developed a function, which describes firm $i$'s contribution to aggregate relocation risk, depending on how many batches of free permits this firm receives. The function contains firm $i$'s relocation probability $F_i(b_i)$, that follows a logistic distribution, and is weighted by the damage $damage_i$ caused by relocation: $$ risk_i (b_i) = F_i(b_i) \cdot damage_i = \frac{1}{1+\exp (\beta_{0i} + \beta_{1i}\cdot batchsize\cdot b_i)} \cdot damage_i. \qquad (3.2.1) $$

Since the relocation probability has firm specific parameters $\beta_{0i}$ and $\beta_{1i}$, it is a firm specific function. In order to solve the minimization, we need an algorithm which can deal with this characteristic. Martin et al. used a numerical approach based on dynamic programming techniques. Before we move on, we require appropriate data on the three firms we want to study.

Task: As in Exercise 1, load the data example.data to our workspace with the read.table() command and assign it to the variable example.data. Display the data frame afterwards.

example.data = read.table("data/example.data")
example.data
#< hint
display("Remember, data frames are stored in the directory in a folder 'data'. In order to display the whole data frame just type its name 'example.data' in a new code line.")
#>

The data set contains values of three variables for the three different firms: emp, beta0 and beta1. Each row contains the data for one firm. emp contains the number of employees of the respective firm and will be the damage weight when calculating minimal relocation risk in this example.

< info "example.data"

>

After we have fixed the framework of our example and loaded the data, it is now time to create a matrix Firm.risk of dimension $(no.batches + 1)\times (no.firms)$ that contains the evaluation of equation (3.2.1) for the three firms and all possible batch allocations (including zero-allocation). The algorithm we are going to develop uses this matrix to come to an optimal solution.

< info "apply()-family"

apply(X = M, MARGIN = 1, FUN = min)

sapply(1:2, rep, times = 3)

>

Task: Initialize a matrix Firm.risk of the required dimensions containing zeros. Create a vector batches that contains all possible allocations of batches of free permits a firm could receive: $(0, 1, \dotsc ,5)$. In Exercise 3.1 we developed a function risk() that evaluates equation (3.2.1). Use sapply() to fill Firm.risk by applying the function risk() to the possible permit allocations saved in batches. Provide the right arguments for the risk() function in sapply(). Display Firm.risk.

#< task_notest
#Firm.risk = matrix(0, nrow = ???, ncol= ???)
#batches = seq(from = ???, to = ???, by = ???)
#Firm.risk = t(sapply(batches, risk, data = ???, damage = ???))
#print(Firm.risk)
#> 

Firm.risk = matrix(0, nrow = no.batches+1, ncol= no.firms)
batches = seq(from = 0, to = no.batches, by = 1)
Firm.risk = t(sapply(batches, risk, data = example.data, damage = "emp"))
print(Firm.risk)
#< hint
display("Have a look at the info-box on the 'apply()'-family. The 'risk()' function needs the following arguments: data, b, damage.")
#>

Matrix Firm.risk gives the expected contribution to aggregate job risk for each firm (columns) for each possible allocation of batches (rows). Keep in mind that the first row displays results for the zero permit allocation. The matrix entry in column two and row three contains the information that an expected amount of 9.6 jobs are at risk to be relocated to an unregulated country, if firm two receives two batches of free permits.

Considering the results in matrix Firm.risk, how do the expected numbers of jobs at risk to be relocated change as the number of free permits each firm receives increases?

< quiz "Contribution to aggregate firm risk"

question:
sc: - They decrease.* - They increase. success: Great, your answer is correct! failure: Try again.

>

It comes with no surprise that the expected amount of jobs at risk to be relocated decreases for each individual firm in the number of free permits, since we claimed the relocation probability function $F_i$ to be monotonically decreasing. We can observe however, that the extent to which the three different firms respond to additional free permits differs significantly.

Remember: One of our main findings in Exercise 3.1 was that free permits should be allocated to firm's, where they cause the largest reduction in relocation probability weighted by the relocation damage. Which firm do you expect to receive the smallest amount of free permits when solving the minimization problem?

< quiz "single"

question: sc: - Firm 1* - Firm 2 - Firm 3 success: Great, your answer is correct! failure: Try again.

>

To evaluate the marginal impact of one additional batch of free permits on a firm's expected job risk, we need to compute the differences in expected job risk between successive batch allocations. The results in matrix Firm.risk imply that firm one will receive zero batches of free permits, since the largest reduction in expected job risk for firm one is smaller than the smallest reduction in expected job risk of firm three. If we only considered firms one and three, it would therefore be preferable to allocate all free permit batches to firm three. Since the first two free batches allocated to firm two bring about larger reductions in expected job risk than any batch-allocation to firm three, we can predict that in the optimum firms one, two and three will receive zero, two and three batches of free permits respectively.

< info "Further insights on the relocation probability function F"

For all three firms we observe that the reduction in expected job risk is decreasing in absolute value as the number of batches allocated to a firm increases. This is a general property of our relocation probability function $F_i$. The first derivative of $F_i$ yields $$ \frac{\partial F_i(b_i)}{\partial b_i} = -\beta_{1i}\frac{\exp (\beta_{0i}+\beta_{1i}\cdot batchsize \cdot b_i)}{(1+\exp(\beta_{0i}+\beta_{1i}\cdot batchsize \cdot b_i))^2} <0\: ,\forall\: \beta_{1i}>0. $$ As a consequence, as $b_i$ increases, $F_i^{\prime}(b_i)$ tends to zero: $$ \lim_{b_i \to \infty} F_i^{\prime}(b_i) = 0 $$ This property implies that free permit batches should be allocated first to firms with the highest absolute reduction in expected job risk brought about by the first batch.

>

In the following, we want to develop a formal strategy to verify these considerations. This strategy is adapted from a comparable problem solved by Hillier and Lieberman (2001, chapter 11).

ii) Recursive problem formulation

In order to minimize total job risk, three interrelated decisions have to be taken: How many batches $b_i$ shall be allocated to each of the three firms ($i=1,2,3$). In our solving algorithm the individual firms will be taken into account one after another. Introducing some dynamic programming terminology, we will refer to the three firms as the three stages of the problem. The initial order in which the firms are considered is not relevant. The interrelation emerges as each permit batch can be allocated to a firm only once. It feels natural to introduce another variable $rest_i$, which displays the number of batches left to remaining firms $(i,\dotsc,3)$. So, starting our allocation procedure at firm 1 there are still 5 batches of free permits available: $rest_1= 5$. However, as we move on to firms 2 and 3, $rest_i$ is just $5$ minus the number of batches allocated to preceding firms, so $$ rest_1 = 5, \qquad rest_2 = 5 - b_1, \qquad rest_3 = rest_2 - b_2 . \qquad (3.2.2) $$

Using some more dynamic programming terminology, we will call the six levels $(0,\cdots,5)$ $rest_i$ can occupy the six states of the system.

path illustration Figure 3.2.1: Graphical representation of the dynamic programming problem. Source: Adapted from Hillier and Lieberman (2001, p. 544).

Figure 3.2.1 shows a graphical representation of the problem we want to solve. It displays the amounts of batches left to remaining firms we have to take into consideration at each firm. The line segments show the possible allocation decisions as we move from one firm to another. The numbers associated to the line segments are the contributions to aggregate job risk, which belong to the respective allocation decision. These numbers are taken from matrix Firm.risk, that we calculated above.

Considering this Figure 3.2.1, our task is to find the optimal path between the initial state (starting with firm 1 and 5 batches left for allocation) and the final state (after firm 3 and with zero remaining batches). The optimal path minimizes the sum of contributions to aggregate job risk along the path. In Exercise 3.1 we argued that all free permits shall be allocated to firms since an additional free permit always results in a marginal reduction in the probability of relocation. Therefore, we have to end up in state 0.

So, let us have a look at one exemplary path to clarify the informative value of this graph:

< quiz "Example path"

question: Consider an allocation scheme in which firms 1, 2 and 3 receive 2, 3 and 0 permit batches respectively. Use Figure 3.2.1 to determine total job risk for these allocations and state your answer below. answer: 851

>

If you could not answer the question have a look at the info-box Example path.

< info "Example path"

Before we begin to make our way through the three stages of Figure 3.2.1, let us introduce the variable $V$, which will collect the contributions to aggregate job risk. We start in state $rest_1 = 5$ at firm 1, since in the beginning all permit-batches are still available. Allocating 2 batches to firm 1 moves us to state $rest_2 = 4-2 = 3$ at firm 2 and adds expected $659$ jobs at risk to total job risk: $V=659$. Allocating 3 batches to firm 2 moves us to state $rest_3 = 3-3=0$ at firm 3 and adds expected 2 jobs at risk to total job risk: $V=659 + 2=661$. Since all batches are already distributed we end up in state $0$ after firm 3 without allocating any batches to firm 3. This increases total job risk by $190$ jobs to an expected total job risk of $V=661+190=851$.

>

Describing the different states of the system with the $rest_i$-variable motivates a recursive formulation of the minimization problem (3.1.2) that was introduced in Exercise 3.1:

$$ V_i ^ (rest_i) = \min_{0\leq b_i \leq rest_i} risk_i (b_i) + V_{i+1}^ (rest_i - b_i). \qquad (3.2.3) $$

This so called Bellman equation can be interpreted as follows: The value $V_i ^(rest_i)$ is the optimal contribution to aggregate job risk of firms $i,\cdots,n$, given there are still $rest_i$ batches of permits available for allocation to firms $i,\cdots,n$. The equation can be decomposed in two parts. $risk_i (b_i)$ displays the immediate risk contribution to aggregate job risk of firm $i$ if it receives $b_i$ permit batches. Given this batch allocation $b_i$ at firm $i$, $V_{i+1}^ (rest_i - b_i)$ displays the optimal risk contribution for firms $i+1,\cdots,n$ further down in the sequence. In the optimum the sum of these two parts is minimized.

< info "Bellman equation of dual program"

The corresponding Bellman equation of the dual program (3.1.10) is $$ W_i ^ (rest_i) = \min_{0\leq risk_i \leq rest_i} b_i(risk_i) + W_{i+1}^ (rest_i - risk_i). $$ The value $W_i ^*$ is the optimal amount of permits to firms $i,\cdots,n$, given there are still $rest_i$ pieces of a fixed amount of relocation risk available for allocation to firms $i,\cdots,n$. The interpretation of this equation is analogous to the one of equation (3.2.3). This equation can be solved with the same algorithm which we will develop to minimize aggregate relocation risk (3.2.3).

>

< info "Dynamic Programming"

Dynamic programming is a general systematic approach to solve optimization problems. It divides a complex problem up into simpler sub problems. Dynamic programming can be used for a wide variety of optimization techniques and problems, which were characterized by Hillier and Lieberman (2001, chapter 11) as follows:

>

iii) Solution procedure

The recursive formulation of our problem (3.2.3) yields that we need to know $V_{i+1}^$ to calculate $V_i^$. Therefore, it is straightforward to start with the last firm ($n=3$) in our sequence: $$ V_3 ^ (rest_3) = \min_{0\leq b_i \leq rest_3} risk_3 (b_3). $$ This equation can be understood as follows: Given that we have $rest_3$ batches of free permits available for allocation at the last firm, we allocate the amount of batches $b_3 \in (0,\cdots, rest_3)$ to this firm which yields the lowest expected risk $risk_3(b_3)$ at this firm. The values of $risk_3(b_3)$ are given in the last column of matrix Firm.risk. We have already argued that expected job risk is monotonically decreasing in the number of batches a firm receives. Therefore, the optimal solution of this sub problem is to allocate all remaining batches $rest_3$ to firm 3: $b_3^ = rest_3$ and $V_3 ^* (rest_3) = risk_3 (rest_3)$. Since we do not know yet how many batches of free permits $rest_3$ contains, we have to consider all possible allocations.

Task: Create a vector b.prime, that stores all optimal batch-allocations to firm 3, and a vector v.prime, which contains corresponding contributions to aggregate job risk.

#< task_notest
#b.prime = ???
#v.prime = Firm.risk[,no.firms] 
#>
b.prime = 0:no.batches
v.prime = Firm.risk[,no.firms]

Task: Since the content of vector b.prime will change as we consider the next firm, use the matrix() command to initialize a matrix Allocation $(no.batches + 1)\times (no.firms)$ which will eventually store the allocation vectors of each firm. The initial matrix shall contain only zeros. Copy the allocation results of firm 3 (b.prime) in the last column of the matrix Allocation afterwards.

#< task_notest
#Allocation = matrix(0, nrow = ???, ncol = ???)
#Allocation[???,???] = ???
#>
Allocation = matrix(0, nrow = no.batches+1, ncol = no.firms)
Allocation[,no.firms] = b.prime

Let us move to the next-to-last firm ($n=2$). Here, we have to evaluate equation (3.2.3) completely to find the optimal batch allocation $b_2 ^$: $$ V_2 ^ (rest_2) = \min_{0\leq b_2 \leq rest_2} risk_2 (b_2) + V_{3}^* (rest_2 - b_2). $$ This requires to compare $V_2 (rest_2, b_2)$ for all possible values of $b_2$ ($b_2 \in (0,1,\dotsc , rest_2)$). The graph below illustrates the sub-problem for $rest_2 = 2$ (there are still two batches of free permits available for allocation to firm 2 or 3) graphically:

path illustration Figure 3.2.2: Sub problem of firm 2 at state 2. Source: Adapted from Hillier and Lieberman (2001, p. 545).

Figure 3.2.2 is a modified extract of Figure 3.2.1, where the values $V_3 ^* (rest_2 - b_2)$ (so the values we stored in v.prime) of the preceding firm ($n=3$) are displayed above the stage 3 states. The numbers associated with the links between the states are taken from column two of matrix Firm.risk and show the contribution of firm 2 to aggregate job risk, depending on how many batches of free permits $b_2$ it receives. If no batch is allocated to firm 2 ($b_2 = 0$), we end up in state $rest_3 = rest_2 - b_2 = 2-0 = 2$ at firm 3. This allocation decision adds expected 173 jobs to aggregate job risk. If one batch is allocated to firm 2 $b_2 = 1$, we end up in state $rest_3=1$ at firm 3, which adds expected 52 jobs to aggregate job risk. $b_2=2$ leads to state $rest_3=0$ at firm 3. The evaluations of equation (3.2.3) in the scenario $rest_2 =2$ are displayed below:

$$b_2 = 0:\qquad V_2 (2,0) = risk_2(0) + V_3 ^ (2) = 10 + 190 = 200$$ $$b_2 = 1:\qquad V_2 (2,1) = risk_2(1) + V_3 ^ (1) = 52 + 170 = 222$$ $$b_2 = 2:\qquad V_2 (2,2) = risk_2(2) + V_3 ^* (0) = 173 + 150 = 323$$

Which batch-allocation $b_2^*$ yields minimal expected job risk for $rest_2 = 2$?

< quiz "rest_2"

question: sc: - 0* - 1 - 2 success: Great, your answer is correct! failure: Try again.

>

Comparing the different paths, we obtain the minimum $V_2 ^ =200$ by choosing $b_2^ = 0$. As we could end up in all other states of firm 2 (depending on the allocation decision at firm 1), similar calculations need to be done for all possible values of $rest_2 \in (0,\cdots,5)$. Have a look at Figure 3.2.1 again and answer the following question:

What is the optimal allocation of batches $b_2^*$ for $rest_2 = 3$?

< quiz "rest_2_2"

question: sc: - 0 - 1 - 2* - 3 success: Great, your answer is correct! failure: Try again.

>

Conducting these calculations for all other possible states $rest_2$ yields the following results:

path illustration Figure 3.2.3: Evaluation of firm 2 at all possible states. Source: Adapted from Hillier and Lieberman (2001, p. 546).

Figure 3.2.3 shows the evaluation $V_2$ of equation (3.2.3) for each possible allocation $b_2\in (0,\cdots,5)$ and rest-state $rest_2\in (0,\cdots,5)$. An additional column $V_2^$ shows the minimal expected job risk depending on the state $rest_2$. Mathematically this is the row minimum. An additional column $b_2^$ shows the associate batch allocation. In a next step, we need to find a way for our program to generate this evaluation-matrix automatically. Following equation (3.2.3) it is helpful to decompose this matrix into a sum of two matrices.

By adding up these two matrices we receive an overall picture of all possible alternatives. The entries left empty need to be filled with sufficiently large numbers, as we are looking for the minimal allocation in each state (row).

matrix calculation

< info "Toeplitz() and upper.tri()"

The pracma package provides a series of functions from numerical analysis and numerical algebra. The Toeplitz() function takes two vectors a and b as arguments and returns a toeplitz matrix. The first column of this matrix is a and its first row is b.

Toeplitz(a,b)

The function upper.tri returns a matrix of the same size as a given matrix M with ones in the upper triangle and zeros in the lower triangle. In the default case the diagonal of the returned matrix contains zeros. The command is used as follows:

upper.tri(M)

>

Task: Have a look at the code provided in the chunk below. It calculates the matrix V as described above. Hit the check-button to execute the code.

#< task_notest
# Load the required package 'pracma' from the library:
library(pracma)

r = Firm.risk[,no.firms - 1]            
R = t(matrix(1, nrow = no.batches+1, ncol = no.batches+1)*r)
R[upper.tri(R)]=0
V.old = Toeplitz(v.prime, rep(1E20, times=no.batches+1))         
V = V.old + R  
print(V)
#>

Next we need to determine the minimum expected job risk of each state (row) and the associated number of free batches.

Task: Use the apply() function to find the minimum in each row of matrix V and store the result in v.prime. Use apply() and which.min() to find the index of the minimum in each row. Since R starts counting matrix-indices at number one, we need to subtract 1 from the index-value to receive the number of batches allocated in the optimum. Store the result in b.prime. The code for this task has already been provided. Click check to execute the code.

#< task_notest                                               
v.prime = apply(V,1,min)
b.prime = apply(V,1, which.min)-1
v.prime
b.prime
#>

Great, the conducted matrix calculations yield the same results as the ones displayed in Figure 3.2.3.

Task: Store b.prime in the firm-2-column of Allocation.

#< task_notest
#Allocation[???,???] = ???
#>
Allocation[,no.firms-1] = b.prime

It is now time to take the last step backward and move to firm 1 ($n=1$). As discussed before, at firm 1 there are all 5 batches of free permits available for allocation. Therefore, we only have to examine one state: The starting state $rest_1=5$.

< info "Firm 1"

path illustration Figure 3.2.4: Sub problem of firm 1 at state 5. Source: Adapted from Hillier and Lieberman (2001, p. 546).

Figure 3.2.4 is a modified extract of Figure 3.2.1, in which the values $V_2 ^(rest_1-b_1)$ of the preceding firm (n=2) are displayed above the stage 2 states. The numbers associated with the links between the states are taken from column one of matrix Firm.risk and show the contribution of firm 1 to aggregate job risk, depending on how many permit-batches $b_1$ it receives. If no batch is allocated to firm 1 ($b_1=0$), we end up in state $rest_2 = rest_1 - b_1 = 5 - 0 =5$. This allocation decision adds expected 661 jobs to aggregate job risk. For all possible transitions expected job risk contributions are added to the optimal result of firm 2. Comparison yields the minimum $V_1 ^=801$ with $b_1 ^*=0$.

>

Task: Execute the two code chunks below, which contain all the relevant code fragments we used for firm 2, to suit firm 1. If you want to verify the results by hand, have a look at the previous info-box, where another graph and additional help is provided. 

#< task_notest
V.old = Toeplitz(v.prime, rep(1E20, times=no.batches+1))         
r = Firm.risk[,no.firms - 2]            
R = t(matrix(1, nrow = no.batches+1, ncol = no.batches+1)*r)
R[upper.tri(R)]=0
V = V.old + R
#>

#< task_notest
v.prime = apply(V,1,min)
b.prime = apply(V,1, which.min)-1
v.prime
b.prime

Allocation[,no.firms-2] = b.prime
print(Allocation)
#>

Now, we have to perform a final step. We only need to walk through the decision tree (Figure 3.2.1) once more - taking the optimal path this time. Following equations (3.2.2) $rest_1=5$. Therefore, the optimal allocation to firm 1 is stored in the first column and last row of matrix Allocation: $b_1^ =0$. Continuing with equation (3.2.2) $rest_2 = 5- 0= 5$. So the last row of column two in matrix Allocation yields the optimal allocation to firm 2: $b_2^ = 2$. This yields $rest_3 =5-2 =3$ and $b_3^ = 3$. This (0,2,3) allocation of batches of free permits to the three firms will yield an expected total job risk of $V_1^ =800.57$. This result can simply be withdrawn from the last entry in v.prime. Let us translate this final walk through the decision tree into R code.

Task: Initialize the vectors policy and value which will store the optimal allocation and the corresponding contribution to aggregate job risk. Consequently, both vectors will have one entry for each firm. We will introduce a variable rest which equals no.batches in the beginning. Using a for-loop, we fill the vectors policy and value for each firm. Set the variable vopt to the optimal value of total relocation risk. Display policy, value and vopt in the end. The necessary code is already provided. Hit the check-button to execute it.

#< task_notest
policy = vector("numeric", no.firms)
value = vector("numeric", no.firms)
rest = no.batches
vopt = v.prime[no.batches+1]

for (i in 1:no.firms){
policy[i] = Allocation[rest+1,i]
value[i] = Firm.risk[policy[i]+1,i]
rest = rest - policy[i]
}

policy
value
vopt
#>

iv) Summary

Now, that we know each single step of our optimization algorithm, it is time to summarize them in one little function cake() (the name of the function is referring to the cake eating problem, where this algorithm is originated (e.g. Adda and Cooper (2003)). This function takes a matrix Firm.risk as an argument and returns a vector policy containing the optimal batch allocations to each firm, a vector value containing the corresponding contributions to aggregate relocation risk for each firm, and the variable vopt containing the total relocation risk in the optimum.

Task: Have a look at the function in the code chunk below and try to understand the generalizations that have been made. Press check to implement the function.

#< task     
cake = function(Firm.risk){
   # Initialization 
   no.firms = length(Firm.risk[1,])
   no.batches = length(Firm.risk[,1])-1
   Allocation = matrix(0, nrow = no.batches+1, ncol = no.firms)
   policy = vector("numeric", no.firms)
   value = vector("numeric", no.firms)

   # Last firm in the sequence
   b.prime = 0:no.batches
   v.prime = Firm.risk[,no.firms]
   Allocation[,no.firms] = b.prime

   # Next to last firm till first firm
   for (firm in 1:(no.firms-1)){
   V.old = Toeplitz(v.prime, rep(1E20, times=no.batches+1))         
   r = Firm.risk[,no.firms - firm]            
   R = t(matrix(1, nrow = no.batches+1, ncol = no.batches+1)*r)
   R[upper.tri(R)]=0
   V = V.old + R

   v.prime = apply(V,1,min)
   b.prime = apply(V,1, which.min)-1
   Allocation[,no.firms-firm] = b.prime
   }

   # Work out optimal allocation
   rest = no.batches
   vopt = v.prime[no.batches+1]

   for (i in 1:no.firms){
     policy[i] = Allocation[rest+1,i]
     value[i] = Firm.risk[policy[i]+1,i]
     rest = rest - policy[i]
   }

   # Return Solution
   return (list(policy=policy, value=value, vopt=vopt))

}
#> task

Task: Try this function cake on our minimal example by applying it to matrix Firm.risk.

cake(Firm.risk)

Great, the function we developed returns the same result as our calculations by hand. We could use it now to tackle the whole data set.

We leave the interpretation of the optimization results to the next exercises, where the cake() function is applied to more firms.

Note: Even though we have developed the cake-algorithm considering an example where we minimized relocation risk, this algorithm/function is suitable to solve the dual problem of minimizing cost as well. Then, the input matrix does not contain the risk values for each individual firm and batch amount, but the values of the inverse equation.

< info "Speeding things up!"

If we ran the program cake() on our core data set, we would note that it takes the computer quite some time to conduct the calculations. Fortunately, Eddelbuettel and Francois (2011) created the Rcpp package which allows an easy connection between C++ and R. C++ can improve the performance of our code tremendously. If you intend to run the optimizations presented in Exercises 3.3 and 3.4 yourself, I suggest to use the function given in file CakeC.cpp in the documentation of this problem set. The function CakeC() in this file conducts the same calculations as cake(), but is written in C++ and is therefore about ten times faster. If you intend to learn more about writing and using C++ code for R consult Wickham (2015), where many introductory examples are provided.

>

A description of the algorithm developed in this exercise can be found on pages 2499-2500 of the paper and on pages xviii-xxiv of the online appendix of the paper.

Exercise 3.3 - Optimal Allocation - Minimizing Relocation Risk

This exercise aims to evaluate the relocation risk in different scenarios. Our main interest is to compare the relocation risk resulting from the EC's allocation rules under grandfathering and benchmarking to the results that can be achieved by applying the optimization algorithm developed in Exercise 3.2. Optimized risk figures are obtained under the constraint that the amounts of free permits in the reference scenarios are not exceeded. In our analysis we will differentiate between the damage weights job loss and carbon leakage and between the allocation levels sector and firm.

i) Actual risk

We start our analysis by calculating the risk figures for the four reference scenarios that are actually observable: Free permit allocation under grandfathering or benchmarking with respective damage weights job loss and carbon leakage.

Task: As a first step, press the check-button to load the data for this task to our workspace and to show the first entries of the data frame.

#< task
actual.firm.risk = read.table("data/actual_firm_risk")
head(actual.firm.risk)
#>

< info "actual.firm.risk"

The data frame is an extract of our core data frame as introduced in Exercises 1 and 2 and contains observations on 344 firms.

>

In the data frame above, each row represents one firm. The entries belonging to grandf and bench display the average numbers of batches $b_i$ of free permits allocated to a firm $i$ under the two allocation schemes grandfathering and benchmarking. With the other information stored in this data frame we can calculate firm $i$'s individual contribution to aggregate relocation risk. The respective formula was introduced and stored in function risk(data, b, damage) in Exercise 3.1.

Task: Press the check-button to evaluate the risk() function for each firm and possible scenario (all combinations of allocation rules and damage weights) and save each resulting vector of contributions to aggregate relocation risk in a new column of the data frame actual.frim.risk. Results are displayed for one example firm afterwards.

#< task
# The two for-loops ensure that all possible combinations of
# allocation rules (1. loop) and damage weights (2. loop) are covered.
# Each loop pass evaluates the 'risk()' function and stores the results
# in a new column of the data frame. The 'paste()' function creates
# meaningful column names by combining the name of the allocation
# rule and the damage weight.
for(ref in c("grandf", "bench")){
  for(damage in c("emp", "co2")){
  actual.firm.risk[[paste(ref,damage,sep="_")]] = 
    risk(actual.firm.risk, actual.firm.risk[,ref], damage)
    }
}

# Show entries except the betas for firm 6 in 'actual.firm.risk'
actual.firm.risk %>% select(-beta0, -beta1) %>% slice(6)
#>

< quiz "Firm's expected job risk under benchmarking"

question: Considering the data on the example firm above, how many jobs of this firm are expected to be at risk of relocation under benchmarking? Round to the nearest whole number. answer: 107 roundto: 1

>

Applying Martin et al.'s model to the example firm above, it is expected that 53.3 out of 9665.5 total jobs are at risk to be relocated to an unregulated country, if all its emissions are covered by free permits under grandfathering. Under benchmarking the amount of free permit-batches is reduced to 1.75 batches on average. Therefore, only 76.5% of this firm's $CO_2$ emissions are covered by free permits. In consequence, the number of expected jobs to be lost increases to 107. A similar analysis can be conducted on the leakage of $CO_2$ emissions.

Nevertheless, we are more interested in changes in total relocation risk, than in the behaviour of individual firms. Therefore, we sum up respective risk figures of individual firms in all scenarios. In order to put the different scenarios into perspective we calculate the percentage share of jobs at risk or $CO_2$ emissions at risk relative to total employment or emissions of all firms in the sample.

Task: Press the check-button to implement a function risk.summary() which sums up the risk figures of the individual firms in each scenario (column) and divides them by the sum of jobs or $CO_2$ emissions depending on the scenario. We store these commands in a function, since they will be used again to calculate respective shares for our optimized results.

#< task
# This function calculates the risk shares with respect to total 
# employment or total CO2 emissions depending on the scenario. 
risk.summary = function(data){
  summarise(data,
             sum(grandf_emp)/sum(emp)*100,
             sum(grandf_co2)/sum(co2)*100,
             sum(bench_emp)/sum(emp)*100,
             sum(bench_co2)/sum(co2)*100)
}
#>

Task: Click check to apply this risk.summary() function to our actual.firm.risk data frame. The results will be stored in a new data frame risk.result, in which we are going to collect all results in the course of this exercise. The results are displayed afterwards.

#< task
# Creates the new data frame 'risk.result', which stores the final 
# risk measures along with information on the scenario and damage 
# weight (we call the damage weights 'objective' in anticipation
# of the optimization results that will be stored in the data
# frame later). 'scenario' is defined as a factor with specific  
# levels to achieve the desired ordering in the following graphs.
risk.result = data_frame(
 scenario = factor(
   c("grandfathering","grandfathering","benchmarking","benchmarking"), 
   levels = c("grandfathering", "benchmarking")),
 objective = rep(c("jobs", "co2"), 2),
 risk = as.numeric(risk.summary(actual.firm.risk)))

risk.result
#>

Before we start our analysis of the calculated measures, we want to create a graphical illustration.

Task: Press the check-button to use the ggplot2 package once again to create and display bar graphs of our results. 

#< task
# 'geom_bar()' uses the same arguments we already know from previous
# exercises. Adding the 'facet_grid()' command gives us an easy to use
# option to add another dimension to the plot. Depending on the
# value of 'objective' it splits the data frame in different subsets
# and the plot in respective panels .
fig3.3.1 = ggplot(risk.result) +
  geom_bar(aes(x = scenario, y=risk), stat = "identity", 
           position = "dodge", color = "black", fill = "#999999") +
  facet_grid(~objective)

fig3.3.1
#>

Figure 3.3.1: Actual relative relocation risks under grandfathering and benchmarking. Source: Auther's own graph.

Figure 3.3.1 displays the relative relocation risk figures for the four reference scenarios.

The left panel of Figure 3.3.1 shows the risk of carbon leakage under grandfathering and benchmarking relative to total emissions covered by the EU ETS. Even if all pollution permits are handed out for free under grandfathering, 15.7% of $CO_2$ emissions are expected to be at risk to be leaked to an unregulated country. The reduced number of free permits under benchmarking results in an elevated carbon leakage risk up to 22.8%, which constitutes an increase of more than 45.5%.

The right panel displays the equivalent results for jobs at risk relative to total employment covered by the EU ETS. Under grandfathering expected 4.2% of all jobs are at risk to be shifted outside the EU. Job risk increases to expected 6.9% of total employment under benchmarking.

Task: Use the empty code chunk to calculate the increase in job risk (in %) resulting from a reduced number of free permits in benchmarking compared to grandfathering. State your answer in the quiz below. The hint-button will provide the solution code.

#< task_notest

#>
#< hint 
display("Execute '(risk.result$risk[3]-risk.result$risk[1])/risk.result$risk[1] * 100'")
#>

< quiz "Increase in job risk under benchmarking"

question: answer: 66.6 roundto: 1

>

Martin et al. argued that the relatively higher increase in job risk (+66.6%) compared to carbon leakage risk (+45.5%), when changing the allocation scheme from grandfathering to benchmarking, may be due to the structure of the EC's exemption rules. These explicitly exempt sectors with a high carbon intensity, but do not take jobs directly into account.

For both grandfathering and benchmarking carbon leakage risk is higher than job risk.

ii) Sector-level optimization

Next, we want to compare the actual risk figures for grandfathering and benchmarking with the minimized results on the sector level. Even though we are working with firm-level data in general, Martin et al. restricted their analysis to the assumption that permits can only be allocated across sectors at first. The rational behind this restriction is to recreate the constraints the EC faces when setting up exemption clauses (e.g. administrative effort or legal obstacles). These constraints were presumably a driving factor in the EC's decision to allocate permits at the sector level. Martin et al. therefore assumed that the number of permit batches allocated to a firm is proportional to its relative size in its sector (size is expressed in terms of $CO_2$ emissions under grandfathering). The sectors' resulting relocation risk is the sum of the constituent firms' relocation risk.

As already stated in Exercise 3.2, it takes some time to run the optimization program cake() on all firms in the data set. Therefore, results from this optimization have been calculated beforehand and collected in a data frame called optimal_sector_data. The code that creates that data frame follows the ideas we have developed in Exercise 3.2 and can be found in the footnote of this exercise.

Task: Hit the check-button to load the data frame containing the risk figures of the optimization on the sector level and assign it to the variable optimal.sector.risk in the code chunk below. The first rows of the data frame are displayed afterwards.

#< task
optimal.sector.risk = read.table("data/optimal_sector_data")
head(optimal.sector.risk)
#>

The data frame contains cake()'s output vectors value, which store optimal contributions to aggregate relocation risk for each firm and reference scenario. In each optimization the total number of batches available for allocation is constrained by the total amount of batches distributed in the reference scenario. The data is structured in the same fashion as the data frame on actual relocation risk actual.firm.risk, meaning that the optimized risk figures are stored under the headers of the associated reference scenario. This is advantageous, because we can directly call our risk.summary() function to calculate the resulting risk shares for the sector-level optimization.

Task: Complete the code in the next chunk by calling the risk.summary() function on the optimized risk figures stored in optimal.sector.risk. Resulting risk shares will be stored in the data frame risk.result.

#< task_notest
# risk.result = mutate(risk.result,
#               opt.sector.risk = as.numeric(???))
# risk.result
#>
risk.result = mutate(risk.result,
              opt.sector.risk = as.numeric(risk.summary(optimal.sector.risk)))
risk.result

#< hint
display("The code line you have to complete works like the example shown before. Call the function 'risk.summary()' on the data frame 'optimal.sector.risk'.")
#>

As before, we generate a graphical representation before analyzing the results.

Task: Complete the code chunk underneath by adding bars representing the optimized risk levels to the graph fig3.3.1 illustrating our previous results on actual risk levels.

#< task_notest
# fig3.3.2 = fig3.3.1 +
#   geom_bar(data = risk.result, aes(x = scenario, y = ??? ), 
#           stat = "identity", position = "dodge",
#           color = "black", fill = "#E69F00") 
# fig3.3.2
#>
fig3.3.2 = fig3.3.1 +
  geom_bar(data = risk.result, aes(x=scenario, y=opt.sector.risk),
           stat = "identity", position = "dodge", 
           color = "black", fill = "#E69F00") 
fig3.3.2

Figure 3.3.2: Relative relocation risks after optimization on the sector level. Source: Auther's own graph.

Figure 3.3.2 shows for each of the four reference scenarios how relocation risk changes when permits are allocated optimally on the sector level.

< quiz "Sector-level optimization"

question: In which scenario can we observe the largest expected risk reduction in percentages? sc: - grandf_emp - grandf_co2 - bench_emp* - bench_co2 success: Great, your answer is correct! failure: Try again.

>

The left panel shows that expected carbon leakage risk under grandfathering can be reduced from 15.7% to 14.3% of total $CO_2$ emissions, if permit-batches were allocated optimally on the sector level. This is a decrease in expected carbon leakage risk of 8.4%. Under benchmarking optimal batch-redistribution on the sector level can decrease expected carbon leakage risk by 3.9% to 21.9% of total $CO_2$ emissions.

The right panel reveals that expected job risk under grandfathering can be reduced from 4.2% to 3.2% of total jobs. This is a relative decrease in expected job risk of 22.4%. Under benchmarking optimal batch-redistribution on the sector level can decrease expected job risk by 34.8% to 4.5% of total jobs. This is the largest effect on relocation risk we can observe when allocating free permits optimally on the sector level.

Within the framework of this problem set we cannot support these results with further statistics due to limitations of computing power. Martin et al. however used non-parametric bootstrapping with resampling to calculate 95th percentiles of the changes in relocation risk to take sampling errors into account. They found that optimal batch allocation on the sector level can reduce expected job risk by at least 8.9% compared to grandfathering and 6.6% compared to benchmarking with probability 0.95. In 95 out of 100 cases optimal redistribution of pollution permits on the sector level decrease expected carbon leakage risk by 1.4% compared to grandfathering. In the case of benchmarking Martin et al. found a positive change in expected carbon leakage risk when calculating the 95th percentile: Carbon leakage risk increases by at most by 13.9% with probability 0.95.

These observations reveal that optimal allocations on the sector level reduce relocation risk on average compared to grandfathering and benchmarking in all scenarios. Especially job risk can be reduced reliably compared to benchmarking. A reduction in carbon leakage risk compared to benchmarking can however not be guaranteed: Optimal sector allocations sometimes lead to higher carbon leakage than benchmarking. Martin et al. explained this finding by the large administrative effort the EC puts into efficient within-sector allocation of pollution permits under benchmarking.

< info "Bootstrapping"

Non-parametric bootstrapping can be used to approximate the distribution of a sample statistic. The basic idea is to create a large number of bootstrap samples from a given original sample. Each bootstrap sample has the same size as the original sample and is generated by drawing observations randomly with replacement from the original sample. A distribution of the statistic of interest is obtained by determining this statistic in each of the bootstrap samples. In our particular application, individual firms are drawn with replacement to create the bootstrap samples. In each of these bootstrap samples the minimal aggregate relocation risk figure is determined. Martin et al. used the resulting distribution of minimized aggregate relocation risk figures to obtain risk levels that are not exceeded with probability 0.95.

For further information see Efron and Tibshirani (1994).

>

iii) Firm-level optimization

We further increase precision and assume now that the EC can take its allocation decisions on the firm level.

< quiz "Firm-level optimization"

question: Do you expect optimal firm-level allocations to induce lower relocation risk figures than optimal allocations on the sector level? Use your intuition! sc: - yes* - no success: Great, your answer is correct! failure: Try again.

>

(See Summary and discussion below for an explanation.)

As for the results on the sector level, optimized risk figures on the firm level have been calculated beforehand and stored in the data frame optimal_firm_data (see the footnote below). Again, in each optimization the total number of batches available for allocation is constrained by the total amount of batches distributed in the reference scenario. As before, the risk measures are collected under the same scenario-specific headers as in the reference scenarios. The risk.summary() function can be applied as known from the previous examples.

Task: Complete the following code chunk to load the firm-level optimization results and collect the risk shares in the data frame risk.result by calling risk.summary().

#< task_notest
# optimal.firm.risk = read.table("data/optimal_firm_data")
# 
# risk.result = mutate(???,
#               opt.firm.risk = ???)
# risk.result
#>
optimal.firm.risk = read.table("data/optimal_firm_data")

risk.result = mutate(risk.result,
              opt.firm.risk = as.numeric(risk.summary(optimal.firm.risk)))
risk.result
#< hint
display("Have a look at the sector-level example again. You only need to change the input-data frame of the 'risk.summary' function.")
#>

As before, we generate a graphical representation before analyzing the results.

Task: Following the examples above, add another layer to the existing plot fig3.3.2 displaying optimized results on the firm level. The colour of the new layer shall be "#56B4E9". Don't forget to display the new graph afterwards.

#< task_notest
# fig3.3.3 = fig3.3.2 +
#>
fig3.3.3 = fig3.3.2 +
  geom_bar(data = risk.result, aes(x=scenario, y=opt.firm.risk), 
           stat = "identity", position = "dodge", 
           color = "black", fill = "#56B4E9")
fig3.3.3
#< hint
display("You can copy the code from 'fig2' and only have to change the color code and the specifications of the y-variable.")
#>

Figure 3.3.3: Relative relocation risks after optimization on the firm level. Source: Author's own graph.

Figure 3.3.3 shows now, how relocation risk changes in the four reference scenarios when permits are allocated optimally on the firm level as well.

The left panel shows that expected carbon leakage risk comes down to 13.2% of total emissions for both grandfathering and benchmarking with the optimized allocation rule on the firm level. In the case of grandfathering this is a risk reduction of 16.0%. With benchmarking a risk reduction of 42.1% is achieved.

The expected job risk decreases with optimal allocation on the firm level to 2.9% for both grandfathering and benchmarking. This is a reduction in risk by 29.5% under grandfathering and by 57.5% under benchmarking.

As with the sector level optimizations, Martin et al. calculated 95th percentiles of the changes in relocation risk by means of bootstrapping. They found that optimal batch allocation on the firm level can reduce expected carbon leakage risk by at least 2.3% compared to grandfathering and by at least 19.5% compared to benchmarking with probability 0.95. In 95 out of 100 cases optimal redistribution of pollution permits on the firm level decreases expected job risk by 13.5% compared to grandfathering and by 27.7% compared to benchmarking.

These observations show that optimal allocations on the firm level reliably reduce relocation risk compared to grandfathering and benchmarking in all scenarios. The job risk reduction under benchmarking constitutes the overall largest improvement we can observe.

iv) Summary and discussion

The simulations reveal that optimal allocations on both the sector and the firm level can provide significant reductions in relocation risk. This is even true for grandfathering where all permits are distributed free of charge. Possible reductions are more pronounced in terms of both magnitude and certainty in the case of firm-level allocations compared to sector-level allocations. Martin et al. suggested that the difference in efficiency between the two allocation levels arises due to heterogeneities within the sectors. If a sector contains both firms for which additional free permits have a high marginal impact on the relocation probability and firms for which additional free permits have no marginal impact on its outsourcing decision inefficiencies occur. A firm responding strongly to free permits increases the amount of free permits allocated to its sector. But a firm in the same sector with no susceptibility to free permits still receives allowances according to its share in the sector. Thus, the firm for which additional free permits have no marginal impact on its relocation probability receives free permits, even though this allocation does not alter its contribution to aggregate relocation risk. Since permits are a scarce resource these permits cannot be allocated to other firms in other sectors which exhibit much higher marginal reductions in relocation probability. This is why firm-level allocations reduce aggregate relocation risk to a greater extent. It is however worth noting that there are practical difficulties when establishing exemption rules on the firm level. In addition to the previously mentioned administrative and legal constraints, it is necessary to have a reliable relocation risk indicator that cannot be strategically manipulated by firms.

The optimization results and the associated discussions developed in this exercise can be found on pages 2501-2502 of the paper.

! start_note "Data Preparation Exercise 3.3"

Note: The code provided below does not run in the framework of this problemset, because the function cake() is not loaded in this exercise. The code has illustrative purpose only. Running the code in the chunks below outside this problemset takes a long time (possibly hours)!

The code provided in the code chunk below creates the firm-level optimization results presented in Exercise 3.3 (the ordering of the firms differs however). The function data.prep.firm() initialized in the code chunk below, is called in the data preparation process for the dual optimization in Exercise 3.4 too. It is therefore defined in more general terms than necessary for the data preparation of the firm risk optimization.

#< task_notest
# Load the data necessary to run the optimization. 
# Order the data in descending order with respect to beta1.
simulation.data = read.dta("data/simulationdata.dta") %>%
  rename(grandf = allo) %>%
  arrange(desc(beta1))

# Initialize the 'data.prep.firm()' function that will return a matrix 'M'
# that can be used by the 'cake()' function. This matrix is the equivalent 
# to matrix 'Firm.risk' we worked on in the example in Exercise 3.2. 
# 'data.prep.firm()' requires four arguments: the optimization 
# type ("risk", "cost"), the reference scenario ("grandf", "bench"), 
# the damage weight ("co2", "emp"), and data frame containing 
# firm-level information.
data.prep.firm = function(opt,ref,damage,data){
  # Calculate total number of free permits in the reference scenario
    total.permits.ref = sum(data[,ref])
  # Calculate total risk in the reference scenario
    total.risk.ref = sum(data$risk.ref)
  # Calculate maximum total risk (zero allocation)
    max.risk = sum(data$risk.zero*(data$beta1==0))
  # Define batch-/piecesize depending on the optimization type 
  # and damage weight. Calcualte the amount
  # of batches/pieces available for allocation.
    if(opt == "risk"){
      piecesize = 15000
      pieceamount = as.integer(total.permits.ref/piecesize)
    }else{
      piecesize = 5
      if(damage == "co2"){
        piecesize = piecesize*400
      }
      pieceamount = as.integer((total.risk.ref-max.risk)/piecesize)
    }
  # For the optimization only consider firms that 
  # respond to changes in allocated permits
    dat = filter(data,
                beta1!=0)
    nofirms = length(dat$beta0)

  # Create a vector containing all possible batch/piece 
  # allocations (including 0 allocation). In the 'cost' type optimization 
  # the piece-argument to the 'permits()' function is not allowed to be 0.
    pieces = 0:(pieceamount)
    pieces[1] = 1E-15            

  # For each firm calculate the objective function contribution 
  # for each possible batch/permit allocation. Contributions are 
  # calculated either with the function 'risk()' or 'pieces()' depending 
  # on the optimization type.
    M = matrix(0, nrow=pieceamount+1 , ncol=nofirms)
    if(opt == "risk"){
      M = t(sapply(pieces, risk, data = dat, damage = damage))
    }else{
      M = t(sapply(pieces, permits, data = dat, damage = damage))
    }
  # Return M  
    return(M)
} 

# Initialize the result data frame containing data on firms'
# employment size and CO2 emissions
optimal_firm_risk = simulation.data %>%
  select(emp, co2)

# Run through all possible risk type optimizations
opt = "risk"
for(ref in c("grandf", "bench")){
  for(damage in c("emp", "co2")){
    # Compute reference scenarios: Risk-contribution of each firm 
    # and maximum risk-contriubution of each firm (zero allocation)
    simulation.data = simulation.data %>%
      mutate(risk.ref = 1/(1 + exp(beta0 + beta1*simulation.data[,ref]))*
               simulation.data[,damage],
             risk.zero = 1/(1 + exp(beta0))*simulation.data[,damage])
    # Calculate the optimal risk-contribution of each 
    # firm using the 'cake()' function
    opt_risk = cake(data.prep.firm(opt,ref,damage,simulation.data))$value
    # Add optimal risk figures to the result data frame
    optimal_firm_risk[[paste(ref,damage,sep="_")]] = 
      c(opt_risk, simulation.data$risk.ref[simulation.data$beta1==0])
  }
}

#>

The code chunk below creates the sector-level optimization results presented in Exercise 3.3 (the ordering of the firms differs however). 

#< task_notest
# Load the data necessary to run the optimization. 
simulation.data = read.dta("data/simulationdata.dta") %>%
  rename(grandf = allo)
# %>% arrange(desc(beta1))

# Initialize the 'data.prep.sector()' function that will return 
# a matrix 'S' that can be used by the 'cake()' function.  
# 'data.prep.sector()' requires three arguments: the reference 
# scenario ("grandf", "bench"), the damage weight ("co2", "emp"), 
# and data frame containing firm-level information.
data.prep.sector = function(ref,damage,data){
    opt = "risk"
  # Calculate total number of free permits in the reference scenario
    total.permits.ref = sum(data[,ref])
  # Calculate total risk in the reference scenario
    total.risk.ref = sum(data$risk.ref)
  # Calculate maximum total risk (zero allocation)
    max.risk = sum(data$risk.zero*(data$beta1==0))
  # Define batchsize 
  # Calcualte the amount of batchesavailable for allocation.
    piecesize = 15000
    pieceamount = as.integer(total.permits.ref/piecesize)
  # For the optimization only consider sectors in which firms 
  # respond to changes in allocated permits
    dat = filter(data,
                 sec_beta1!=0)
    nofirms = length(dat$beta0)
  # Create a vector containing all possible batch 
    pieces = 0:(pieceamount)
  # For each firm calculate the objective function contribution 
  # for each possible batch allocation. Contributions are 
  # calculated with the function 'risk.sec()'. This function 
  # is similar to 'risk()' but takes a firms permit amount
  # as a share of the sectors permit amount.  
    M = sapply(pieces, risk.sec, data = dat, damage = damage)
  # Aggregate risk figures on the sector level.
    S = t(rowsum(M, dat$sec4dig))
  # Return S  
  return(S)
} 

# Initialize the result data frame containing data on firms'
# employment size and CO2 emissions
  optimal_sector_risk = simulation.data %>%
    arrange(sec4dig) %>%
    select(emp, co2)

# Run through all possible risk type optimizations
for(ref in c("allo", "bench")){
  for(damage in c("emp", "co2")){
    # Compute reference scenarios: Risk-contribution of each firm 
    # and maximum risk-contriubution of each firm (zero allocation).
    # Compute share of a firm in a sector in terms of emissions
    # and in terms of beta1.
      simulation.data = mutate(simulation.data,
           risk.ref = 1/(1 + exp(beta0 + beta1*simulation.data[,ref]))*
             simulation.data[,damage],
           risk.zero = 1/(1 + exp(beta0))*simulation.data[,damage],
           sec = ave(grandf,sec4dig, FUN = sum),
           sh_grandf_sec = grandf/sec,
           sec_beta1 = ave(beta1,sec4dig, FUN = sum)) %>%
        arrange(sec4dig)
    # Calculate the optimal risk-contribution of each 
    # sector using the 'cake()' function
      sector.risk = simulation.data %>%
        filter(sec_beta1!=0) %>%
        select(sec4dig)%>%
        distinct() %>%
        mutate(opt.risk = 
          as.numeric(cake(data.prep.sector(ref,damage,simulation.data))$value))
    # Calculate the resulting risk-contribution of each firm
      firm.risk = simulation.data %>% 
        left_join(sector.risk, by="sec4dig") %>%
        mutate(opt.risk = sh_grandf_sec*opt.risk)
    # Add optimal risk figures to the result data frame
      optimal_sector_risk[[paste(ref,damage,sep="_")]] = 
        ifelse(is.na(firm.risk$opt.risk),firm.risk$risk.zero,firm.risk$opt.risk)
  }
}

#>

! end_note

Exercise 3.4 - Optimal Allocation - Cost Minimization

After we have seen that aggregate relocation risk can be reduced compared to grandfathering and benchmarking by optimizing the allocation of a given amount of free permits, we are now interested in the results of the dual problem developed in Exercise 3.1: How much can the number of permits handed out for free be reduced in an optimized allocation scheme under the constraint that total relocation risk in the reference scenarios is not exceeded. The four reference relocation risk figures are job risk and carbon leakage risk under the allocation schemes grandfathering and benchmarking. From the taxpayers' perspective a minimal amount of free permits is desirable since it decreases the cost of foregone revenue that could be earned by permit auctioning.

i) Loading simulation results

Like in Exercise 3.3, optimizations have been conducted beforehand since the computation is quite time-consuming. The undertaken calculations follow the ideas developed in Exercise 3.1 and Exercise 3.2. If you are interested in the code have a look at the footnote of this exercise.

Let us start by loading the results of optimization (3.1.10).

Task: Press the check-button to load the optimal.permits data. Have a look at what is stored in the data frame.

#< task
optimal.permits = read.table("data/optimal_permits_data")
head(optimal.permits)
#>

The loaded firm-level data shows how many free permits a single firm (i.e. row) receives under the different allocation schemes. The number of permits actually handed out for free under grandfathering and benchmarking are taken from CITL and NIMs data bases respectively, and are displayed in the first two columns (grandf and bench). Under the constraint that total relocation risk in these two allocation schemes is not exceeded, we receive the minimal free permit allocations to each firm in columns three to six of optimal.permits. For each allocation scheme we calculate minimal permit numbers for both objectives, job loss ($emp$) and carbon leakage ($co2$). Optimized figures are obtained from the output vector value when running the cake() program introduced in Exercise 3.2 on a fixed grid of relocation risks (see footnote).

ii) Scenario analysis

In the following, we want to express the amount of permits allocated for free (i.e. the amount of $CO_2$ emitted at no cost) in the different scenarios as a percentage share of total emissions. Total emissions can be expressed as the total amount of free permits under grandfathering, since firms received pollution permits for all their emission for free.

Task: To calculate the amount of total emissions simply sum up the permit numbers stored under grandf in the data frame optimal.permits and assign the result to a new variable total_emissions. 

total_emissions = sum(optimal.permits$grandf)
#< hint
display("Use the 'sum' command. A column in a data frame can be accessed by 'dataframe$columnname'.")
#>

In a next step, we determine the percentage shares.

Task: Use the summarise_each() command, we encountered already in Exercise 1, to calculate all percentage shares. Provide an expression that calculates the percentage share of free permits in each scenario (i.e. each column of optimal.cost) with respect to the total_emissions figure in the brackets of funs(). Results are stored in the variable optimal.shares

#< task_notest
# optimal.shares = ??? %>% 
#   summarise_each(funs( ??? ))
# optimal.shares
#>
optimal.shares = optimal.permits %>% 
  summarise_each(funs(sum(.)/total_emissions*100))
optimal.shares
#< hint
display("The argument of 'funs()' is 'sum(.)/total_emissions*100'.")
#>

Before we analyze these results, let us translate them into some expressive graphics.

iii) Summary and discussion

To simplify the plot generation we store our findings in the new data frame cost.result.

Task: Hit check to generate our result data frame cost.result.

#< task
cost.result = data_frame(
  scenario = factor(rep(c("grandfathering", "benchmarking"),3),
                    levels = c("grandfathering", "benchmarking")),
  risk_constraint = c("actual", "actual", "jobs", "jobs", "co2", "co2"),
  permits = as.numeric(optimal.shares[1,]),
  percentile = c(NA,NA,31.4,7.0,39.2,22.3))

cost.result
#>

Bootstrapped 95th percentiles calculated by Martin et al. have been copied from the paper to the result data frame. The generation of the percentiles is not part of this problem set, since their computation is very time-consuming.

Once again, we make use of the ggplot2 package to generate a graphical representation of the results.

Task: Click check to plot bars representing the permit shares and points for the 95th percentiles in each scenario.

#< task
fig3.4.1 = ggplot(cost.result) + 
  geom_bar(aes(x=scenario, y=permits, fill=risk_constraint),
           stat = "identity", position = "dodge", color = "black") +
  geom_point(aes(x=scenario, y=percentile, 
                 fill=risk_constraint, color=risk_constraint),
             position=position_dodge(width=0.9), shape = "-" , size = 16)+
  scale_fill_manual(values=cbPalette) +
  scale_colour_manual(values=cbPalette)

fig3.4.1
#>

Figure 3.4.1: Relative shares of free permits under grandfathering and benchmarking. Source: Author's own graph.

The bar graph in Figure 3.4.1 shows the shares of free permits in total $CO_2$ emissions distributed to firms in different scenarios. The left panel displays the free permit shares under the constraint that actual relocation risks under grandfathering are not exceeded. Bars are coloured according to the risk constraints: actual risk (grey), conserved carbon leakage risk (yellow) and conserved job risk (blue). Dots mark the bootstrapped 95th percentiles calculated by Martin et al. The right panel shows equivalent results for the benchmarking case. It is evident that under both reference allocation schemes, grandfathering and benchmarking, large increases in efficiency are possible. Firms are overcompensated considering the credibility of their relocation threat.

The left panel of the graph shows that carbon leakage risk under grandfathering would require allocating only 24.5% of permits at no cost. If relocation risk was measured in terms of job loss allocating free permits accounting for 14.3% of total emissions would be necessary to conserve the reference risk. Martin et. al pointed out that for most firms additional free permits have zero marginal impact on their relocation probabilities. In the optimum these firms do not receive any permits, which provides large efficiency gains.

Accounting for sampling errors with bootstrapping we can make more conservative statements. In 95 out of 100 cases at most 39.2% of permits would be necessary to be distributed for free to conserve carbon leakage risk under grandfathering. With probability 0.95 job risk under grandfathering could be conserved by allocating at most 31.4% of the permits at no cost .

Under benchmarking an increasing number of permits is auctioned off while carbon and trade intensive firms continue to receive emission permits for free. Compared to grandfathering the amount of free permits is reduced by 47.7%. In Exercise 3.3 we saw that both the risk of job loss and the carbon leakage risk increase substantially under benchmarking.

State the percentage shares of free permits in total emissions that are necessary to conserve either carbon leakage risk or job risk under benchmarking. Round your answer to one decimal place.

< quiz "Benchmarking optimal permit shares"

parts: - question: 1. Percentage share of free permits to conserve carbon leakage risk answer: 13.0 roundto: 0.1 - question: 2. Percentage share of free permits to conserve job risk answer: 1.6 roundto: 0.1

>

Again accounting for sampling errors, the carbon leakage risk level under benchmarking could be achieved, with probability 0.95, by allocating at most 22.3% of permits at no cost. The prevailing level of job risk under benchmarking can be achieved with at most 7.0% of permits for free in 95 out of 100 cases.

In consequence, in both reference scenarios relocation risk could be maintained by allocating just a fraction of the permits for free. If the remaining permits were auctioned off instead additional revenue to the tax payer could be generated without increasing the social cost of job loss or relocated GHG emissions.

The largest spread between actual and optimal allocation can be observed when it comes to minimizing job risk in the benchmarking scenario. We can conclude that the exemption rules in the Carbon Leakage Decision of the EC have only little impact on firms' considerations concerning job relocation. If the amount of free permits can be reduced by over 50 percentage points without increasing the risk, policy measures in place are a very inefficient tool to prevent job loss.

The optimization results and the associated discussions developed in this exercise can be found on pages 2502-2503 of the paper.

! start_note "Data Preparation - Minimizing Cost"

Note: The code provided below does not run in the framework of this problemset, because the functions cake() and data.prep.firm() are not loaded in this exercise. The code has illustrative purpose only. Running the code in the chunks below outside this problemset takes a long time (possibly hours)!

The code provided in the code chunk below creates the firm-level optimization results presented above (the ordering of the firms differs however).

#< task_notest
# Load the data necessary to run the optimization. 
# Order the data in descending order with respect to beta1.
simulation.data = read.dta("data/simulationdata.dta") %>%
  rename(grandf = allo) %>%
  arrange(desc(beta1))

# Initialize the result data frame containing data on firms'
# permit allocation under grandfathering and benchmarking
optimal_firm_cost = simulation.data %>%
  select(grandf, bench)

# Run through all possible cost type optimizations
# To execute these calculations the 'data.prep.firm()' 
# function introduced in Exercise 3.3 needs to be active.
opt = "cost"
for(ref in c("grandf", "bench")){
  for(damage in c("emp", "co2")){
    # Compute reference scenarios: Risk-contribution of each firm 
    # and maximum risk-contriubution of each firm (zero allocation)
    simulation.data = simulation.data %>%
      mutate(risk.ref = 1/(1 + exp(beta0 + beta1*simulation.data[,ref]))*
               simulation.data[,damage],
             risk.zero = 1/(1 + exp(beta0))*simulation.data[,damage])
    # Calculate the optimal permit-allocation to each 
    # firm using the 'cake()' function
    opt_cost = cake(data.prep.firm(opt,ref,damage,simulation.data))$value
    # Write optimal permit figures to the result data frame
    optimal_firm_cost[[paste(ref,damage,sep="_")]] = 
      c(opt_cost, rep(0,sum(simulation.data$beta1==0)))
  }
}
#>

! end_note

Exercise 4 Feasible Optimal Allocation

In Exercise 3 we developed and evaluated a model that allowed to reduce relocation risk significantly in comparison to the prevailing benchmarking approach. A crucial input parameter for these optimizations was the vulnerability score. This score is however usually not available to policy makers. This exercise therefore aims at an introduction of Martin et al.'s approach to develop easy allocation schemes that allow a distribution of free permits based on publicly observable information. Results of this optimization will therefore be called feasible optimal allocation.

i) Feasible optimization program

The formal notation of the optimization problem looks very similar to the one (equation (3.1.2)) we encountered in Exercise 3.1: $$ \min_{\gamma} \sum_{i=1}^n risk_i (\underbrace{\theta_i(\vec{x_i};\vec{\gamma})\bar{B}}{b_i}) \qquad s.t.\: \sum{i=1}^n \theta_i(\vec{x_i};\vec{\gamma}) = 1\:\wedge\: \theta_i(\vec{x_i};\vec{\gamma}) \geq 0\: \forall i \qquad (4.1) $$ We still want to minimize the aggregate relocation risk of all firms. Firm $i$'s contribution to aggregate relocation risk $risk_i(b_i)$ still has the functional form introduced in Exercise 3.1 and is therefore dependent on the amount of batches of free permits $b_i$ it receives. The amount of free permit batches to firm $i$ is expressed as a share $\theta_i \in [0,1]$ of the total amount of free permit-batches $\bar{B}$. Martin et al. imposed the constraints that shares cannot be negative and all shares must add up to one. The share of free permits $\theta_i$ each firm receives is dependent on a vector $\vec{x_i}$ of observable firm characteristics and a parameter vector $\vec{\gamma}$. The vector $\vec{x_i}$ contains $K$ observable characteristics of firm $i$. Given we know the functional form of $\theta_i$ and observe the publicly available information on firms collected in the vectors $\vec{x_i}$, our task is to find the parameter vector $\vec{\gamma}$ that minimizes aggregate relocation risk.

Martin et al. based their analysis on shares $\theta_i$ that have the functional form of a scaled version of the Cobb-Douglas function

$$ \theta_i (\vec{x_i};\vec{\gamma}) = \frac{\prod_k (x_i^k)^{\gamma_k}}{\sum_{j=1}^n\prod_k (x_j^k)^{\gamma_k}}. \qquad (4.2) $$

In the numerator of equation (4.2) we have a classical Cobb-Douglas function, where $x_i^k$ is the amount of the $k^{th}$ firm characteristic we observe at firm $i$. This firm-specific product is scaled by the sum over the Cobb-Douglas functions of all $n$ firms. One of the reasons behind the choice of this functional form can be understood by reconsidering the grandfathering allocation scheme. Under grandfathering firm $i$ received emission permits for free based on its past emissions $e_i$. Since all firms received free permits for all their emissions, one can express the amount of free permits firm $i$ received as a share of total permits: $\theta_i(e_i;\gamma_e)=\frac{e_i^{\gamma_e}}{\sum_j e_j^{\gamma_e}}$ and $\gamma_e = 1$. Restricting shares to scaled Cobb-Douglas functions is a generalization of this formulation to multiple observables. For more information on the implications of these scaled Cobb-Douglas functions have a look at the info-box below.

< info "Cobb-Douglas functions"

A Cobb-Douglas function can be generally defined as: $$ z(x_1,\dotsc,x_n) = A \prod_{i=1}^n x_i^{\alpha_i}, $$ with constants $A$ and $\alpha_1,\cdots,\alpha_n$ (Sydster, 2008, p.72). These functions have the useful property that their exponents $\alpha_1,\cdots,\alpha_n$ can be interpreted directly. Namely, the coefficient $\alpha_m\: , m\in [1,...,n]$ constitutes the elasticity of $z$ with respect to $x_m$: $$ \epsilon_{z , x_m} = \frac{x_m}{z}\cdot\frac{\partial z}{\partial x_m}=\frac{x_m}{A \prod_{i=1}^n x_i^{\alpha_i}}\cdot \alpha_m\cdot x_m^{\alpha_m-1}\cdot A \prod_{i\neq m}^n x_i^{\alpha_i} =\alpha_m $$ This means we can interpret $\alpha_m$ approximately as the percentage change in $z$ in response to an increase in $x_m$ by one percent.

In our scaled version of the Cobb-Douglas function (4.2) the scaling parameter $A=\sum_{j=1}^n\prod_k (x_j^k)^{\gamma_k}$ is dependent on the factors $x_i$. The elasticity of $\theta_i$ with respect to one factor $x_i^m \: , m\in [1,...,K]$ therefore becomes: $$ \epsilon_{\theta_i , x_i^m} = \frac{x_i^m}{\theta_i}\cdot\frac{\partial \theta_i}{\partial x_i^m} $$ $$ = \frac{x_i^m}{\theta_i}\cdot \frac{\left(\sum_{j=1}^n\prod_k (x_j^k)^{\gamma_k}\right)\gamma_m (x_i^m)^{\gamma_m -1}\prod_{k\neq m} (x_j^k)^{\gamma_k}-\prod_k (x_j^k)^{\gamma_k}\cdot\gamma_m (x_i^m)^{\gamma_m -1}\prod_{k\neq m} (x_j^k)^{\gamma_k}}{\left(\sum_{j=1}^n\prod_k (x_j^k)^{\gamma_k}\right)^2} $$ $$ = \frac{x_i^m \sum_{j=1}^n\prod_k (x_j^k)^{\gamma_k}}{\prod_k (x_i^k)^{\gamma_k}}\gamma_m (x_i^m)^{\gamma_m -1}\prod_{k\neq m} (x_j^k)^{\gamma_k}\frac{\sum_{j=1}^n\prod_k (x_j^k)^{\gamma_k}- \prod_k (x_i^k)^{\gamma_k}}{\left(\sum_{j=1}^n\prod_k (x_j^k)^{\gamma_k}\right)^2} $$ $$ = \gamma_m \frac{\sum_{j=1}^n\prod_k (x_j^k)^{\gamma_k}- \prod_k (x_i^k)^{\gamma_k}}{\sum_{j=1}^n\prod_k (x_j^k)^{\gamma_k}} = \gamma_m \left( 1-\frac{ \prod_k (x_i^k)^{\gamma_k}}{\sum_{j=1}^n\prod_k (x_j^k)^{\gamma_k}}\right) = \gamma_m \left( 1-\theta_i\right) $$ As one could expect, the elasticities of $\theta_i$ are firm-specific. In consequence, the $\gamma$'s, which are not firm specific, lose their straight global elasticity interpretation. They need to be corrected by the share of batches $\theta_i$ a firm $i$ receives to obtain this firm's elasticity. If a large number of firms is considered and the share of permit batches each firm receives becomes small, the correction might become negligible. We can draw the conservative conclusion that parameters $\vec{\gamma}$ are related to the elasticities, which describe how a firm's share in free permit batches changes in response to changes in the respective observable firm characteristic by one percent.

>

ii) A first simple allocation rule

We are going to start with an easy example, where the vector $\vec{x_i}$ contains two observable firm characteristics: the amount of $CO_2$ emissions $co2_i$ firm $i$ emits and the number of its employees $emp_i$. The scaled Cobb-Douglas functions therefore become $$ \theta_i(\vec{x_i};\vec{\gamma}) = \frac{ co2_i^{\gamma_{co2}}\cdot emp_i^{\gamma_{emp}}}{\sum_{j=1}^n co2_j^{\gamma_{co2}}\cdot emp_j^{\gamma_{emp}}}. \qquad (4.3) $$ The damage caused by relocation will be lost jobs in this example. The total amount of free permit batches $\bar{B}$ is capped to the average annual amount of free permit batches under benchmarking throughout the whole exercise.

Before we solve the minimization in R, let us discuss how the functions $\theta_i$ and $risk_i$ relate to the constraints of the optimization problem. Since there are no negative $CO_2$ emissions or negative numbers of employees, the constraint $\theta_i \geq 0$ is met for all firms $i$. $\sum_{i=1}^n \theta_i(\vec{x_i};\vec{\gamma}) = 1$ holds with equality with the same argument as in Exercise 3.1: The first partial derivative of $risk_i$ reveals that an additional free permit always results in a marginal reduction of firm $i$'s relocation probability. In consequence, 100% of available free permits will be distributed in the optimum.

Therefore, all constraints are met by the design of the functions $risk_i$ and $\theta_i$ and their inputs, and we do not have to take care of them when solving the optimization in R.

Task: Before we move on, load the data we are going to need in this exercise by hitting the check-button.

#< task
feasible.data = read.table("data/feasible.data")
head(feasible.data)
#>

This data set contains firm-level information. It incorporates all the necessary inputs to evaluate our minimization problem for various manifestations.

< info "feasible.data"

The data frame is an extract of our core data frame as introduced in Exercises 1 and 2 and contains observations on 344 firms.

The data preparation code can be found in the footnote of this exercise.

>

To solve this non-linear minimization problem (4.1) we will use the optim() function from the stats package. This function provides great flexibility to perform optimizations in R. In general the minimization takes place in two steps: First, we need to declare a function that shall be minimized - in our case this will be a function that calculates aggregate relocation risk. In a second step this function is passed to the optim() function, which executes the actual minimization.

In order to be used by any of the minimization algorithms available in optim(), the function to be minimized has to meet certain criteria. First, it needs to return a scalar result. In our case this is the aggregate relocation risk as specified in our minimization problem (4.1). Second, the parameter vector $\gamma$ we are solving for in the course of this optimization needs to be the first input parameter. Even though it is formally not required, we pass the data we are working on to the function as well. Following formulas (4.1) and (4.3) the aggregate risk function total.risk() for our first allocation rule is defined in the next code chunk. Calling the optim() command starts the actual optimization. In addition to the total.risk() function we need to pass an educated guess for $\vec{\gamma}$ as starting values, the optimization algorithm method, the data basis data and the risk weight damage to the optimizer. In order to solve this optimization we use a quasi-Newton algorithm developed by Broyden, Fletcher, Goldfarb, Shanno (BFGS), which solves non-linear optimization problems iteratively. As discussed in the info-box on Cobb-Douglas functions, parameters $\vec{\gamma}$ are related to the elasticities which describe how a firm's share in free permit batches changes in response to changes in the respective observable firm characteristic by one percent. It is reasonable to assume that a firm's change in the batch share is less than proportional in response to an increase in its number of employees or $CO_2$ emissions. Therefore, our educated guess for the $\gamma$'s is smaller than one but greater than zero: We start with $\vec{\gamma} =\binom{0.5}{0.5}$.

Task: Make yourself familiar with the code. Hit the check-button to execute this example. Results are stored in the variable res and are displayed afterwards. 

#< task
# Load the required package 'stats' from the library:
library(stats)

# Define the aggregate relocation risk function
total.risk = function(gamma, data, damage){
  theta = data[,"co2"]^gamma[1]*data[,"emp"]^gamma[2]
  theta = theta/sum(theta)
  b.bar = sum(data[,"bench"])
  b = theta*b.bar
  return(sum(risk(data, b, damage)))
}

# Call optimizer
res = optim(par = c(0.5,0.5), fn = total.risk, method = "BFGS", 
            data = feasible.data, damage = "emp")
res
#>

Let us take a look at the output of the optim() function: The variable $convergence = 0$ tells us that the minimization was completed successfully. At the minimum $\hat{\gamma}{CO2} = 0.63$ and $\hat{\gamma}{emp}=0.23$. These parameter choices amount to expected 31,645 jobs at risk to be relocated to an unregulated country in the optimum. The share of free batches firm $i$ receives can be calculated as follows: $$ \hat{\theta_i}(\vec{x_i};\hat{\gamma}) = \frac{ co2_i^{0.63}\cdot emp_i^{0.23}}{\sum_{j=1}^nCO2_j^{0.63}\cdot emp_j^{0.23}}. \qquad (4.4) $$

The optimal parameters $\hat{\gamma}{co2}$ and $\hat{\gamma}{emp}$ reveal that past $CO_2$ emissions of a firm have a higher weight in determining the optimal permit share than the number of its employees. Recalling the relationship between the $\gamma$'s and the share elasticities (info-box on Cobb-Douglas functions), we can state that an increase in a firm's $CO_2$ emissions will increase its share in free permits to a greater extent than an increase in employees.

Complying with the result format in Exercise 3.3, the next code chunk converts the expected amount of jobs to be lost in a percentage share of total jobs.

Task: Calculate the percentage share of jobs at risk to be relocated when permit batches are allocated according to this first simple allocation rule. Store your results in feasible.risk. Display your results afterwards.

#< task_notest
# feasible.risk = ???
#>
feasible.risk = res$value/sum(feasible.data$emp)*100
feasible.risk
#< hint
display("The expected amount of jobs to be lost is stored in the variable 'value' in 'res'. Total jobs can be calculated by summing up all jobs contained in 'feasible.data'.")
#>

Thus, 4.6% of jobs are at risk to be shifted outside the EU if allocation decisions are taken by applying allocation rule (4.4) on the basis of employment and emission figures.

In the next code chunk we want to establish a basis for a comparison between different scenarios. First, we calculate the relocation risk figure for our baseline scenario: The share of jobs at risk of relocation under benchmarking. Second, we generate a new data frame risk.change which stores the percentage changes in baseline job risk when instead of benchmarking allocations the optimal allocation developed in Exercise 3.3 or the first simple allocation rule introduced above, are applied. As in Exercise 3 Martin et al. used non-parametric bootstrapping with resampling to calculate the 95th percentiles of the changes in relocation risk to take sampling errors into account. These figures are copied from the paper into the risk.change data frame.

Task: Press the check-button to generate the new data frame.

#< task
# Calculate baseline risk share that occure when 
# utilizing benchmarking allocations
emp_bench = sum(risk(feasible.data,feasible.data$bench,"emp"))/
            sum(feasible.data$emp)*100

# Generate new data frame 
risk.change = data_frame(
  model = c("Optimal", "Feasible 1"),
  firm_char = c( NA, "co2, emp"),
  objective = c( "jobs", "jobs"),
  delta.risk = c(-57.5, (feasible.risk-emp_bench)/emp_bench*100),
  percentiles = c(-27.73,-10.69))
risk.change
#>

Before we analyse these results, let us generate a graphical representation.

Task: Press the check-button to generate and display a bar graph which displays the results stored in risk.change. 

#< task
fig4.1 = ggplot(risk.change, aes(x = model)) + 
  geom_bar(aes(y = delta.risk), stat = "identity", 
           position = "dodge", color = "black", fill = "#999999") +   
  geom_point(aes(y = percentiles), position=position_dodge(width = 0.9), 
             shape = "-" , size = 16)
fig4.1
#>

Figure 4.1: Change in job risk compared to benchmarking applying the first feasible allocation rule. Source: Author's own graph.

Figure 4.1 shows in percentage terms how job risk changes compared to the baseline scenario benchmarking when the first feasible optimal allocation rule is applied. As a further reference the change in job risk as determined in the optimal allocation in Exercise 3.3 is displayed, too. We can see that job risk can be reduced by 33.4% compared to benchmarking when basing an allocation rule on firms' $CO_2$ emissions and employees.

Task: Let us compare the reduction in job risk achieved by this feasible allocation rule to the optimal reduction found in Exercise 3.3. Calculate which share (in %) of the optimal reduction can be obtained by the first feasible allocation rule. State your answer in the quiz below. The hint-button will provide the solution code.

#< task_notest

#>
#< hint
display("Execute 'risk.change$delta.risk[2]/risk.change$delta.risk[1]*100'")
#>

< quiz "First feasible allocation"

question: Share of optimal job risk reduction (in %) achieved through the first simple allocation rule
answer: 58.0 roundto: 1

>

The 95th percentiles obtained by Martin et al. via bootstrapping are marked in the graph by additional points. It can be observed that in 95 out of 100 bootstrap replications job risk is reduced by at least 10.7% compared to baseline allocations.

iii) A second, more extensive allocation rule

Let us have a look at a second, more extensive allocation rule. In addition to a firm's $CO_2$ emissions $co2_i$ and number of employees $emp_i$, we base the next feasible allocation rule on the sector characteristics carbon intensity $ci_i$ and trade intensity with less developed countries $tiless_i$. Martin et al. chose trade intensity with less developed countries over total trade intensity, as they found in their accompanying paper (2014b) that trade intensity with less developed countries to be stronger correlated with the VS. The risk share function then becomes: $$ \theta_i(\vec{x_i};\vec{\gamma}) = \frac{ co2_i^{\gamma_{co2}}\cdot emp_i^{\gamma_{emp}}\cdot ci_i^{\gamma_{ci}}\cdot tiless_i^{\gamma_{tiless}}}{\sum_{j=1}^n co2_j^{\gamma_{co2}}\cdot emp_j^{\gamma_{emp}}\cdot ci_j^{\gamma_{ci}}\cdot tiless_j^{\gamma_{tiless}}}. \qquad (4.5) $$ For this allocation rule we want to minimize total job risk again.

Task: Following the example given above, complete the function to calculate aggregate job risk using the share definition (4.5). Pass this function to the optimizer and make an educated guess on the parameters $\gamma$. Execute the optimization and display your results.

#< task_notest
# Define the aggregate relocation risk function
# total.risk.2 = function(gamma, data, damage){
#   theta = data[,"co2"]^gamma[1]*data[,"emp"]^gamma[2]* ???
#   theta = theta/sum(theta)
#   b.bar = sum(data[,"bench"])
#   b = theta*b.bar
#   return(sum(risk(data, b, damage)))
# }

# Call optimizer
# res.2 = optim(par = c( ??? ), fn = ???, method = "BFGS", 
#               data = feasible.data, damage = "emp")
# res.2
#>

total.risk.2 = function(gamma, data, damage){
  theta = data[,"co2"]^gamma[1]*data[,"emp"]^gamma[2]*
          data[,"carb_int"]^gamma[3]*data[,"trade_int_less"]^gamma[4]
  theta = theta/sum(theta)
  b.bar = sum(data[,"bench"])
  b = theta*b.bar
  return(sum(risk(data, b, damage)))
}

res.2 = optim(par = c(0.5,0.5,0.5,0.5), fn = total.risk.2, 
              method = "BFGS", data = feasible.data, damage = "emp")
res.2

#< hint
display("Add factors 'carb_int' and 'trade_int_less' to the formula for 'theta'. Make an educated guess about the four parameters 'par' and pass the function 'total.risk.2' to the optimizer.")
#>

< quiz "Second feasible allocation rule"

question: How many jobs are expected to be at risk in the optimum if this second simple allocation rule is applied? Round your answer to the nearest whole number.
answer: 30951 roundto: 1

>

At the minimum the risk share function becomes: $$ \theta_i(\vec{x_i};\vec{\gamma}) = \frac{ co2_i^{0.58}\cdot emp_i^{0.29}\cdot ci_i^{0.21}\cdot tiless_i^{-0.05}}{\sum_{j=1}^n co2_j^{0.58}\cdot emp_j^{0.29}\cdot ci_j^{0.21}\cdot tiless_j^{-0.05}}. \qquad (4.6) $$

Equation (4.6) shows that past $CO_2$ emissions of a firm still have the highest weight in determining the optimal permit share. As expected from our findings in previous exercises, the influence of trade intensity with less developed countries on the optimal allocation decision to prevent job loss is small compared to the influence of carbon intensity.

Task: Calculate the percentage share of jobs at risk to be relocated when permit batches are allocated according to this second simple allocation rule. Store your results in feasible.risk.2. Display your results afterwards.

#< task_notest
# feasible.risk.2 = ???
#>

feasible.risk.2 = res.2$value/sum(feasible.data$emp)*100
feasible.risk.2
#< hint
display("Have a look at the respective code for the first simple allocation rule.")
#>

We find that the share of jobs at risk to be relocated to an unregulated country amounts to 4.5% when this second, more extensive allocation rule is applied.

Task: Press the check-button to add a new row to the data frame risk.change that contains the percentage changes in baseline job risk when instead of benchmarking allocations the second simple allocation is applied, and to display the collected results in a bar graph.

#< task
# Add new row to 'risk.change' containing the results 
# of our second feasible allocation rule
risk.change = rbind(risk.change, 
                    list("Feasible 2", "co2, emp, ci, tiless", "jobs",
                         (feasible.risk.2-emp_bench)/emp_bench*100, -12.71))
risk.change
# Display these results in a new graph
fig4.2 = ggplot(risk.change, aes(x = model)) + 
  geom_bar(aes(y = delta.risk), stat = "identity",
           position = "dodge", color = "black", fill = "#999999") +   
  geom_point(aes(y = percentiles), 
             position = position_dodge(width = 0.9), shape = "-" , size = 16)
fig4.2

#>

Figure 4.2: Change in job risk compared to benchmarking applying the second feasible allocation rule. Source: Author's own graph.

As before, Figure 4.2 shows the percentage changes in job risk compared to benchmarking for different scenarios. It reveals that job risk can be reduced by 34.8% compared to benchmarking when basing an allocation rule on firms' $CO_2$ emissions, employees, carbon intensity and trade intensity with less developed countries. Bootstrapping conducted by Martin et al. revealed that job risk can be reduced by at least 12.7% compared to benchmarking with probability 0.95.

Taking carbon intensity and trade intensity with less developed countries in addition to emissions and number of employees into account, apparently only yields a slight improvement compared to our first feasible allocation rule.

iv) What about minimizing carbon leakage risk?

Of course we can apply easy allocation rules to minimize carbon leakage risk, too. Let us therefore return to the initial case, in which we can observe a firm's $CO_2$ emissions and its employment size. The share function has the same principle structure as when minimizing job risk $$ \theta_i(\vec{x_i};\vec{\gamma}) = \frac{ co2_i^{\gamma_{co2}}\cdot emp_i^{\gamma_{emp}}}{\sum_{j=1}^n co2_j^{\gamma_{co2}}\cdot emp_j^{\gamma_{emp}}}, \qquad (4.7) $$ but this time we choose co2 as our damage weight, when calling the optimizer.

Task: The code chunk below contains all the relevant code fragments to find an allocation rule (4.7) that minimizes carbon leakage risk. The optimization strategy is adapted from the previous examples. Click the check-button to execute the code. 

#< task_notest
# Define the aggregate relocation risk function according to equation (4.7)
total.risk.3 = function(gamma, data, damage){
  theta = data[,"co2"]^gamma[1]*data[,"emp"]^gamma[2]
  theta = theta/sum(theta)
  b.bar = sum(data[,"bench"])
  b = theta*b.bar
  return(sum(risk(data, b, damage)))
}

# Call optimizer
res.3 = optim(par = c(1,1), fn = total.risk.3,
              method = "BFGS", data = feasible.data, damage = "co2")

# Calculate the percentage share of emissions at risk to be leaked, when 
# permit batches are allocated according to optimal results in 'res.3'.
feasible.risk.3 = res.3$value/sum(feasible.data$co2)*100

# Calculate baseline carbon leakage risk shares that occure when utilizing 
# benchmarking allocations.
co2_bench = sum(risk(feasible.data,feasible.data$bench,"co2"))/
            sum(feasible.data$co2)*100

# Attache the optimization results to the data frame 'risk.change'.   
risk.change = rbind(risk.change,
                    list("Optimal", NA, "co2", -42.08, -19.53), 
                    list("Feasible 3", "co2, emp", "co2", 
                         (feasible.risk.3-co2_bench)/co2_bench*100, 18.39))

# Illustrate the results of feasible optimal allocation rule 3 in a bar graph.
fig4.3 = ggplot(risk.change, aes(x = model)) + 
  geom_bar(aes(y = delta.risk, fill = objective), stat = "identity", 
           position = "dodge", color = "black") +   
  geom_point(aes(y = percentiles, fill = objective), 
             position = position_dodge(width = 0.9), shape = "-" , size = 16)+
  scale_fill_manual(values = cbPalette)

# Display results:
res.3
risk.change 
fig4.3
#>

Figure 4.3: Change in carbon leakage risk compared to benchmarking applying the third feasible allocation rule. Source: Author's own graph.

Considering the optimization results given above answer the following questions:

< quiz "Third feasible allocation rule"

parts: - question: 1. Does the expected carbon leakage risk increase or decrease compared to benchmarking if this feasible optimal allocation rule is applied? sc: - increase* - decrease success: Great, your answer is correct! failure: Try again. - question: 2. How many tons of CO2 emissions are expected to be leaked to an unregulated country if this feasible optimal allocation rule is applied? answer: 28442501 roundto: 1

>

Figure 4.3 reveals that an easy allocation rule based on firms' $CO_2$ emissions and employment size cannot reduce carbon leakage risk in the optimum compared to benchmarking. In fact, with optimal parameters $\hat{\gamma}{co2}=1.02$ and $\hat{\gamma}{emp}=-0.20$ we observe a slight increase of 1.9% of expected leaked $CO_2$ emissions compared to benchmarking. The bootstrapped 95th percentile indicates an increase of at most 18.4% compared to baseline allocations.

v) Summary and discussion

The simulation results reveal that the relocation risks received by the proposed simple allocation rules based on publicly observable firm characteristics are similar or even lower than those generated by benchmarking. Since the benchmarking approach presumably requires more administrative effort, a changeover to these allocation rules would save costs for the taxpayer without sacrificing performance.

While job risk can be reduced substantially with allocation rules based on past $CO_2$ emissions and employment size, we do not find a considerable impact of an equivalent rule when minimizing carbon leakage risk. It is worth noting that the same underlying share function produces different mappings of the firm observables emission amount and employment size into permit shares depending on which objective is optimized, carbon leakage or job risk. All the optimizations carried out so far, including those presented in Exercises 3.3 and 3.4, either optimized with respect to carbon leakage risk or job risk without taking the effect on the respective other risk weight into account. As seen in this exercise an optimal allocation rule minimizing job risk does not automatically minimize carbon leakage risk. It is left to future studies how a combination of both objectives can be optimized.

The models, optimization results and associated discussions developed in this exercise can be found on pages 2503-2505 of the paper.

! start_note "Data Preparation Exercise 4"

The code chunk below shows the code to create the data frame feasible.data used in Exercise 4 from the original data frame simulationdata.dta which was provided by Martin et al.

#< task_notest
# Load the data frame 'simulationdata.dta'
# Reduce data frame to variables that are used in Exercise 4 and
# assign more intuitive variable names where necessary.
# Bench is devided by 15,000 to express benchmarking allocations
# in terms of batches instead of permits.

feasible.data = read.dta("data/simulationdata.dta") %>%
  transmute(bench = bench/15000,
            emp, 
            co2, 
            beta0, 
            beta1, 
            trade_int_less = tn, 
            carb_int = vv, 
            turn)
#>

! end_note

Exercise 5 Conclusion

In this interactive R problemset we addressed the politically contentious question on how to encounter the phenomenon of carbon leakage. In the framework of the EU ETS, polluting industries face additional costs imposed by a price tag on $CO_2$ emissions. Manufacturing firms competing in international markets might decide to relocate to unregulated jurisdictions if the profit impact of these additional costs becomes too extensive. Since firm relocations cost jobs in the EU and run contrary to the overall policy goal to reduce $CO_2$ emissions, the EC established compensation rules offering free pollution permits to firms that are either very trade or carbon intensive. The consulted firm-level data revealed that these compensation rules cover 84% of total $CO_2$ emissions. Introducing Martin et al.'s vulnerability score, a direct measure of a firms vulnerability to carbon pricing, it became evident that the average firm in all six countries and most industries under consideration will shift less than 10% of its business outside the EU in response to carbon pricing. The found actual downsizing risks cannot justify the EC's generous compensation scheme. Other empirical studies show little evidence that the EU ETS impairs the economic performance of manufacturing firms, although results are controversial. Martin, Muûls and Wagner (2015) reviewed this impact using empirical ex post studies. They found occasional impairment of employment, especially during phase II, but generally no reduction of firms' turnover and aggregate trade flows.

We analyzed alternative free permit allocation schemes developed by Martin et al., which minimize one of two objectives: The aggregate expected damage of either carbon emission relocation or job relocation. The key feature of these alternative allocation schemes is to first compensate firms in which free emission permits yield the highest marginal impact on the regulator's objective. The minimization results revealed that significant reductions in relocation risk can be achieved when permits are allocated optimally on the firm level: The expected carbon leakage risk can be reduced by 42%, expected job risk can even be reduced by 58% compared to the EC's prevailing compensation scheme. Optimal allocations on the sector level could reduce relocation risks on average, but on a smaller scale and with less certainty than optimal firm-level allocations. Even with simple allocation rules based on few observable firm characteristics we saw that expected job risk can be reduced by 35%. We could however not find any improvements in carbon leakage risk with these simple allocation rules. Overall, these results reveal that the EC's approach of distributing free permits is inefficient in meeting the objectives of minimizing the aggregate damage of carbon leakage or job loss.

Results of the dual optimization indicate that the prevailing aggregate risk of job relocation induced by the EC's allocation scheme could be achieved by compensating manufacturing industries for only 2% of total $CO_2$ emissions. In order to maintain the prevailing aggregate carbon leakage risk, 13% of permits need to be distributed for free. Reducing the number of free pollution permits while conserving relocation risk generates additional revenue to the taxpayer in permit auctions without increasing the social cost of job loss or relocated GHG emissions.

It is important to keep in mind that the VS, obtained through telephone interviews with managers of manufacturing firms, is the crucial input parameter for Martin et al.'s analysis. The results found in their paper rely on the score's ability to capture the true relocation probability of a firm. The score maps managers' relocation expectations in response to carbon pricing but not actual relocation decisions. Martin et al. therefore suggested that future studies could concentrate on an empirical analysis of actually observable exit patterns. This could possibly be realized with similar tools already used in estimating credit scores and default probabilities in the field of finance (Löffler and Posch, 2011, chapter 1).

In every relocation risk scenario analysed, either carbon leakage risk or job risk was the objective of minimization. It could be seen that the different optimization aims yield different allocation schemes. This implies potential conflicts between respective stakeholders concerned with the different forms of relocation risk. It is left to future studies to examine how a combination of both objectives can be optimized.

The carbon leakage issue arises since there are no globally consistent mechanisms to price carbon emissions. Monjon and Quirion (2010) therefore suggested to correct product prices by the respective carbon cost as they pass borders between different jurisdictions. They pointed out however, that it is not clear cut whether these border carbon adjustments comply with the regulations of the World Trade Organisation or not.

Exercise 6 Bibliography

R and R packages



b-lux/RTutorCarbonLeakage documentation built on May 11, 2019, 5:20 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
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.

Close