# The Carbon Premium: Correlation or Causation? Evidence from S&P 500
Companies

Namasi G. Sankar  Indian Institute of Science Education and Research Pune,
Pashan, Pune - 411008, Maharashtra, India. E-mail:
<EMAIL_ADDRESS>Suryadeepto Nag  Indian Institute
of Science Education and Research Pune, Pashan, Pune - 411008, Maharashtra,
India, and Faculty of Geosciences and Environment, Institute of Geography and
Sustainability, University of Lausanne, Geopolis, Mouline, CH-1015 Lausanne,
Switzerland. E-mail: suryadeepto.nag<EMAIL_ADDRESS>Siddhartha P. Chakrabarty  Department of Mathematics, Indian Institute of
Technology Guwahati, Guwahati-781039, Assam, India. E-mail<EMAIL_ADDRESS>Sankarshan Basu  Department of Finance and Accounting, Indian Institute of
Management Bangalore, Bengaluru-560076, Karnataka, India. E-mail:
<EMAIL_ADDRESS>
###### Abstract

In the context of whether investors are aware of carbon-related risks, it is
often hypothesized that there may be a carbon premium in the value of stocks
of firms, conferring an abnormal excess value to firms’ shares as a form of
compensation to investors for their transition risk exposure through the
ownership of carbon instensive stocks. However, there is little consensus in
the literature regarding the existence of such a premium. Moreover few studies
have examined whether the correlation that is often observed is actually
causal. The pertinent question is whether more polluting firms give higher
returns or do firms with high returns have less incentive to decarbonize? In
this study, we investigate whether firms’ emissions is causally linked to the
presence of a carbon premium in a panel of 141 firms listed in the S&P500
index using fixed-effects analysis, with propensity score weighting to control
for selection bias in which firms increase their emissions. We find that there
is a statistically significant positive carbon premium associated with Scope 1
emissions, while there is no significant premium associated with Scope 2
emissions, implying that risks associated with direct emissions by the firm
are priced, while bought emissions are not.

JEL Classification: G10, G32

Keywords: Carbon premium, Carbon risk, Sustainable finance, Transition risk,
Green finance

##  1 Introduction

A pertinent question in sustainable finance is whether financial markets are
pricing carbon risk. Studies often find that firms with higher carbon
footprints are associated with excess stock returns that cannot be explained
by other relevant factors. Does this imply that investors are demanding higher
returns to stay invested in firms that could have large long-term carbon risk
exposure? Or is it possible that firms that perform well and offer good
returns to investors have less incentive to be environmentally responsible?
There are also a host of studies that find conflicting results about the
existence of said excess returns (Chakrabarty and Nag, 2023). While there are
considerable differences in the empirical specifications used in various
studies, very few studies make an attempt to establish a robust, causal effect
between greenhouse gas (GHG) emissions and stock prices. The primary objective
of this paper is to investigate whether the association between firms’ carbon
footprints and their stock prices is causal, or merely a correlation.

According to a report published by CDP (Griffin, 2017), just 50 companies
contributed to half the GHG emissions in the world in 2015. Therefore, if we
are to not exceed ecological boundaries, firms may need to reduce their energy
consumption from fossil fuels in a substantial manner. Such a transition,
whether it comes from government regulation such as emissions trading (Zhang
and Wei, 2010), carbon taxes (Rutherford, 1992; Hamilton and Cameron, 1994;
Oswald et al., 2023), carbon rationing (Fawcett, 2004; Howell, 2012) or
voluntary initiatives from the firms, will have negative consequences on the
revenue and profit generation, through various channels such as costs of
equity capital (Kim et al., 2015), stranded assets (Curtin et al., 2019), or
through indirect effects of reduced consumption (Semieniuk et al., 2021).
Entire sectors may exhibit declines (Di Cosmo and Hyland, 2013), and even
switching to cleaner production may not achieve positive impacts immediately
(Cucchiella et al., 2017). CDP’s report also states that of the 224 fossil
fuel companies in the sample, 30% were public-investor owned and 11% were
private-investor owned, implying that low-carbon transitions could imply large
risks for investors. Therefore, it is imperative to study whether investors
are aware of such risks and whether this manifests as a risk premium or a
penalty in stock prices.

One possibility is that investors are already aware of risks, and have a
preference for stocks of companies that are environmentally responsible.
Pástor et al. (2021) argue that in equilibrium, these (green) firms would have
lower expected returns than polluting (brown) firms, as the stock prices of
green firms may be inflated due to investor preferences. However, in the
presence of a growing consensus on impending climate change, real markets can
hardly be expected to be in equilibrium, and Pastor et al. (2022) find that
there exists an abnormal surplus realized return in green firms rather than
brown firms, which the authors attribute to an ongoing increase in the demand
for green stocks. In contrast, Bolton and Kacperczyk (2021b) find a positive
association between emissions and stock returns and conclude that markets have
already priced carbon risk to a certain extent. The authors interpret the
association as a premium demanded by investors to continue investing in risky,
carbon-intensive firms. They explain that, since carbon-intensive stocks may
have significant future risks, an investor would invest in it only if the firm
were providing greater short-term gains that may offset the long-term
exposure. This implies that if carbon risk were appropriately priced, there
would be a “fair” carbon premium such that the expected abnormal returns would
be equal to the expected long-term exposure (Nag et al., 2023). Thus, firms
are incentivized to provide an “abnormal” short-term premium to investors to
keep them invested in their shares, which results in a positive carbon
premium.

Despite risks, however, Bolton and Kacperczyk (2021b) find unexceptional
divestments in polluting firms. Krueger et al. (2020) find that institutional
investors, although aware of the financial risks related to climate change, do
not think of divestments as an optimal strategy and tend to opt for risk
management instead. Risk management, unlike divestments, demands careful
analyses and relies more on measures that could adequately quantify the risks.
Thus, the existence of a carbon premium would be of particular interest to
academics and investors interested in developing strategies for hedging
climate risk (Andersson et al., 2016; Roncalli et al., 2020, 2021).

In general, there is a lack of consensus in the empirical literature with
regard to the existence of a positive carbon premium. In a systematic review,
Chakrabarty and Nag (2023) find positive associations between GHG emissions
and stock returns in 9 out of 26 studies surveyed, the remaining finding
neutral, mixed, or even negative results.111Nine articles –Oestreich and
Tsiakas (2015); Del Prado (2019); Wen et al. (2020); Alessi et al. (2021);
Bolton and Kacperczyk (2021a, b); Walkshäusl (2021); Gurvich and Creamer
(2022), and Park et al. (2022) – estimate positive carbon premia or a negative
greenium. Seven articles – Ravina and Hentati Kaffel (2019); Kuang and Liang
(2020); Bernardini et al. (2021); Ji et al. (2021); Muñoz (2021); Meyer (2022)
and Pastor et al. (2022) – estimate positive greeniums or negative premia. Six
articles – Maddox (2018); Zhang et al. (2018); Görgen et al. (2020); Brunborg
and Haldorsen (2021); Aswani et al. (2022); Vymyatnina and Chernykh (2022) and
Zhang (2022) – find no premia. The remaining four studies – Mazzarano et al.
(2021); Santi and Moretti (2021); Witkowski et al. (2021) and Lo et al. (2022)
report mixed results. This lack of consensus in the literature could, in part,
be explained by different policies in different markets studied. For instance,
several studies that survey European markets (Oestreich and Tsiakas, 2015;
Ravina and Hentati Kaffel, 2019; Bernardini et al., 2021) would estimate risks
in the context of the European Union’s Emissions Trading Scheme, where the
presence of carbon allowances may induce effects in the market which may not
be replicable elsewhere, especially in markets such as those of the United
States or other emerging economies that do not have such regulations in place.
Strangely, findings do not agree even when the study is set in the same market
circumstances. Oestreich and Tsiakas (2015) finds a statistically significant
positive effect associated with the dirty-minus-clean portfolio, while
Bernardini et al. (2021) and Ravina and Hentati Kaffel (2019) find the
opposite effect. Chakrabarty and Nag (2023) also find that there may exist
methodological biases in determining whether a study reports a positive carbon
premium or otherwise, and the necessity of more studies that explore causality
in the estimation of carbon premia.

We use publicly available data on GHG emissions and company financials of 141
S&P 500-listed firms traded in US markets to analyze whether a larger volume
of carbon emissions leads to higher investor returns at present. Our final
sample of 141 firms together constitute about a third of the total market
capitalization of US markets. We study stock price movements between 2018 and
2019 so that inferences are not confounded by market shocks induced by the
COVID-19 pandemic in 2020 and 2021. If the value of a stock is thought of as
the demand for the stock, the carbon premium can be thought of as the carbon
elasticity of demand. Studying the elasticity of demand tells us how stock
prices respond to a one percent change in GHG emissions. Studying the demand
curves for the shares of a company against the environmental price would
require knowledge of a company’s potential returns for a variety of GHG
emission levels, which cannot be estimated in a market that is not in
equilibrium, with public sentiment constantly changing. Accordingly, we study
a panel of firms so as to study the demand for shares at a market level. This
poses additional challenges related to confounding variables and firm-fixed
effects. To address these, we employ a difference-in-differences
specification, controlling for a variety of financial parameters that could
independently influence stock prices.

A crucial assumption in difference-in-difference estimation is that of random
treatment (Roth et al., 2023), which is unlikely to be the case in firms’ GHG
emissions. Bolton and Kacperczyk (2021b) find that a company’s carbon
emissions are determined by various factors such as size, return on equity,
etc. Under such circumstances, whether a company increases its emissions or by
how much a company increases its emissions, i.e., the “treatment”222Treatment
here is a change in GHG emissions, while control is no change. Since changes
in firm emissions can be any real number, we essentially have multiple
treatment categories or rather continuous treatment doses. would be affected
by various observables, not controlling for which would result in selection
bias in the assignment to treatment and control groups (Tucker, 2010).
Endogeneity in treatment assignment may make it difficult to make causal
inferences as the effect of GHG emissions may also be a manifestation of the
effects of the factors that affect a firm’s GHG emissions. Very few studies in
the literature have controlled for selection bias in the past. Wen et al.
(2020) and Santi and Moretti (2021) both use matching techniques to address
selection bias. One drawback of matching techniques is that not all treated
firms need to have a corresponding control firm with similar pre-treatment
characteristics. This leads to a contraction of the sample size. Furthermore,
even if pre-treatment covariates are similar, they may vary differently in the
study period. Therefore, following Hirano and Imbens (2004), we weight our
firms in the regressions with the inverse of their propensity scores, to
ensure that the assignment to treatment, given covariates, is ”as good as
random”. Weighting by the inverse of propensity scores simply controls for
selection bias by balancing the representation of pre-treatment covariates in
treatment and control groups, and since a regression follows, differences in
changes of other variables can be controlled for.

To quantify firms’ emissions, we use Scope 1 and Scope 2 emissions. Scope 1
emissions are those emissions that are directly produced by the firm on its
premises or on assets controlled by the firm. Scope 2 emissions are the
emissions from purchased electricity. These emissions are not produced, but
rather consumed by the firm and hence add to the footprint. We do not include
Scope 3 emissions because these are reported less frequently and include
various different types of emissions, the contributions of which to a premium
if there is one, would not necessarily be the same. However, one type of Scope
3 emissions, namely, direct use-phase emissions, is an extremely important
category. These are emissions that arise from the products sold by the firm.
This is particularly relevant for companies that produce, store or transport
fossil fuels. In the event of a low-carbon transition, fossil fuel companies’
revenues would be hurt not only because of the emissions they consume but also
from reduced consumption of the emissions they sell. We incorporate a variable
in our model that proxies for sold emissions of fossil fuel companies, and
estimates the risks associated with sold emissions.

From our analysis, we find a statistically significant, positive carbon
premium associated with Scope 1 emissions between 2018 and 2019, although we
do not find such a premium for Scope 2 emissions. We find that firms with more
assets, fewer liabilities, and low existing emissions were most likely to
increase their Scope 1 emissions. Our estimates after weighting by propensity
scores, in order to remove selection bias, strengthens the estimate of the
carbon premium in both size and significance, implying that the carbon premium
is larger for firms that are less likely to increase their emissions. We also
find that carbon premia are not symmetric, and that while increases in
emissions result in increases in an abnormal increase in stock prices, a
decrease in emissions does not lead to a decrease in stock prices. This may be
because carbon risk is already valued at the market, and the market expects
production to become substantially cleaner, and the pace at which firms are
reducing their emissions may not be fast enough. However, an increase in
emissions may be unexpected among growing concerns on climate change, and
leads to an increase in the valuation of transition risk for said firms,
causing investors to demand greater returns on the short term.

The paper is organized as follows. In the following section, we present the
data sources and the descriptive statistics of the data under consideration.
In Section 3, we present the methodology, with the results and main
conclusions being presented in Section 4 and 5, respectively.

##  2 Data

We use free, publicly available data from S&P 500-listed firms for our
analysis. Based on availability, we collected emissions data for 209 companies
listed in US markets, which are a part of the S&P 500 index. This data was
collected from voluntary company disclosures, either to third-party
organizations such as CDP, which the companies themselves also published or
from ESG reports of the companies. The Greenhouse Gas Protocol classifies
carbon emissions into three categories, namely, Scope 1, Scope 2, and Scope 3.
Scope 1 emissions include emissions produced by the firm on its own assets,
while Scope 2 emissions include emissions from purchased electricity. The
latter are not emissions produced by the company, but rather adding to the
company’s carbon footprint. Scope 3 is a more heterogenous category that
includes a variety of different types of emissions such as emissions from the
consumption of services such as those from business-related travel to the
emissions arising from the products manufactured by companies. We restrict our
analysis to Scope 1 and Scope 2 emissions. This is for two major reasons.
First, the number of firms that have data available on Scope 1 and Scope 2
emissions greatly outnumber the number of firms that have data on Scope 3
emissions (in addition to Scope 1 and Scope 2 emissions) and if we were to
include firms that had data on all three types of emissions, our sample would
shrink considerably. Secondly, while Scope 1 emissions clearly include only
emissions produced by the firm, and Scope 2 emissions include only emissions
consumed by the firm, this is not true for Scope 3 emissions, which could
include both, and there is little reason to believe that any risk would be
homogeneously linked to different kinds of Scope 3 emissions. This makes the
interpretation of any quantitative results difficult, and a proper analysis
would require more clarity while looking within Scope 3 emissions. However,
this is beyond the scope of this work. Most of the companies in our database
had data for Scope 1 and Scope 2 emissions for the years 2018 and 2019, and we
selected firms that had data for both years.

| Means (standard deviation)
---|---
| 2018 |  | 2019 |  | $\Delta$
GICS | Log Scope 1 | Log Scope 2 | Log Total |  | Log Scope 1 | Log Scope 2 | Log Total |  | Log Scope 1 | Log Scope 2 | Log Total
Sector | Emissions | Emissions | Emissions |  | Emissions | Emissions | Emissions |  | Emissions | Emissions | Emissions
All Sectors | 12.6607 | 12.3149 | 13.7497 |  | 12.6715 | 12.1249 | 13.7022 |  | 0.0108 | -0.1900 | -0.0476
| (3.1822) | (2.1451) | (2.4059) |  | (3.1178) | (2.3150) | (2.4012) |  | (0.3198) | (1.0537) | (0.3072)
Energy | 15.8383 | 13.9467 | 16.0241 |  | 15.7922 | 13.6064 | 15.9206 |  | -0.0462 | -0.3403 | -0.1035
| (1.5379) | (1.4383) | (1.4946) |  | (1.5201) | (1.4846) | (1.5067) |  | (0.2428) | (0.6967) | (0.2532)
Materials | 14.6101 | 14.0287 | 15.1028 |  | 14.6046 | 13.9630 | 15.0806 |  | -0.0056 | -0.0657 | -0.0222
| (1.5744) | (1.3040) | (1.4564) |  | (1.5856) | (1.3180) | (1.4688) |  | (0.0635) | (0.0741) | (0.0521)
Industrials | 13.5263 | 12.2232 | 14.2680 |  | 13.4607 | 12.1402 | 14.1862 |  | -0.0656 | -0.0830 | -0.0818
| (2.6229) | (1.4218) | (2.0105) |  | (2.6653) | (1.3666) | (2.0525) |  | (0.1182) | (0.1350) | (0.0987)
Consumer Discretionary | 12.2436 | 13.4256 | 13.7077 |  | 12.2387 | 13.3919 | 13.6832 |  | -0.0049 | -0.0338 | -0.0245
| (1.3523) | (1.1724) | (1.2055) |  | (1.3346) | (1.1955) | (1.2176) |  | (0.0933) | (0.0485) | (0.0434)
Consumer Staples | 13.3777 | 13.4810 | 14.2956 |  | 13.3615 | 13.0902 | 14.0090 |  | -0.0162 | -0.3907 | -0.2866
| (1.8078) | (1.9157) | (1.9099) |  | (1.8213) | (1.3867) | (1.6306) |  | (0.0371) | (0.8564) | (0.6171)
Health Care | 11.2774 | 11.3483 | 12.3596 |  | 11.2949 | 11.3375 | 12.3321 |  | 0.0175 | -0.0109 | -0.0274
| (1.5621) | (2.5807) | (1.3661) |  | (1.4997) | (2.5098) | (1.3814) |  | (0.1542) | (0.2349) | (0.1761)
Financials | 8.5815 | 10.2833 | 10.5032 |  | 8.5374 | 10.1018 | 10.3614 |  | -0.0441 | -0.1815 | -0.1418
| (1.2514) | (0.6718) | (0.7265) |  | (1.2627) | (0.7571) | (0.7653) |  | (0.1737) | (0.2104) | (0.1343)
Information Technology | 9.9745 | 11.5278 | 12.0444 |  | 10.2007 | 11.6759 | 12.1787 |  | 0.2262 | 0.1481 | 0.1343
| (2.9704) | (2.2317) | (1.8110) |  | (2.8408) | (2.1638) | (1.9971) |  | (0.6995) | (0.8332) | (0.5739)
Communication Services | 13.8350 | 15.7170 | 15.8588 |  | 13.8064 | 15.5690 | 15.7274 |  | -0.0286 | -0.1480 | -0.1314
| (0.0000) | (0.0000) | (0.0000) |  | (0.0000) | (0.0000) | (0.0000) |  | (0.0000) | (0.0000) | (0.0000)
Utilities | 15.9515 | 12.3496 | 16.3178 |  | 15.8799 | 11.4556 | 16.2275 |  | -0.0717 | -0.8940 | -0.0903
| (2.2097) | (1.9429) | (1.5875) |  | (2.1826) | (3.2777) | (1.5731) |  | (0.1285) | (2.3904) | (0.0922)
Real Estate | 10.2411 | 12.0469 | 12.3568 |  | 10.2890 | 11.9954 | 12.3250 |  | 0.0479 | -0.0515 | -0.0318
| (0.6985) | (1.7374) | (1.4855) |  | (0.6354) | (1.6382) | (1.3712) |  | (0.0970) | (0.1048) | (0.1198)

Table 1: Sector wise GHG emissions of firms in 2018 and 2019. Includes 141
firms from 11 sectors. Averages are unweighted.

For stock price data, we use data from Yahoo Finance, and for company
financials, we use data from Google Finance, which has data on various
variables of interest such as revenue, expenses, profits, earnings per share,
etc. Both the databases, namely, Yahoo Finance and Google Finance are free and
publicly available. Based on the availability of Scope 1 emissions, Scope 2
emissions, and data on a variety of financial variables, we construct a panel
of 141 firms that maximizes the number of variables that we can use for our
analysis without compromising severely on the size of our sample. To classify
firms by sectors, we use the Global Industry Classification Standard developed
by MSCI and S&P Dow Jones indices, which classifies industries into 11 sectors
– Energy, Materials, Industrials, Consumer Discretionary, Consumer Staples,
Health Care, Financials, Information Technology, Communication Services,
Utilities, and Real Estate.

(a) Distribution of Changes in Log Emissions

(b) Distribution of Changes in Emissions Intensity

Figure 1: The distributions of changes in emissions (log) and changes in
emissions intensity between 2018 and 2019. Emissions intensity is calculated
by dividing GHG emissions by firm revenues of the same year. Includes 141
firms from 11 sectors.

Table 1 presents the Scope 1, Scope 2, and total (Scope 1 + Scope 2) emissions
for firms from all sectors. Overall, emissions were reduced between 2018 and
2019, with an average reduction of 4.65% in aggregate emissions. However, most
of the reduction is from Scope 2 emissions, which declined by close to 18%. On
the contrary, Scope 1 emissions marginally increased by 1.08% over 2018
levels. Within sectors, the biggest aggregate reduction in emissions is in the
Consumer Staples sector with close to a 25% reduction in emissions, while
Materials and Health Care industries observe only marginal reductions of under
3% each. While most sectors also reduced their Scope 1 emissions, the
Information Technology sector increased their Scope 1 emissions by over 25%
above their 2018 level, probably due to increased production, which biases the
average for all sectors combined. An implication of the large fall in Scope 2
emissions, against the small rise in Scope 1 emissions is that companies
reduced their energy consumption substantially. The smaller reduction in Scope
1 emissions implies that the production of emissions from companies changes
little in comparison (to bought emissions). This may be a consequence of
businesses opting for more responsible electricity consumption as a means of
reducing their carbon footprints, as opposed to significantly improving their
production mechanisms. Figure 1 displays the distribution of changes in
aggregate emissions (subfigure 1(a)) and emissions intensities (subfigure
1(b)) of firms in the sample.333Emissions intensity is calculated by scaling
aggregate annual emissions by the revenue generated in that year. In both
cases, we find tall peaks at the center, with few outliers, suggesting that
most firms have comparable emission and electricity consumption patterns.
However, the distribution of emissions intensity is considerably more skewed
with a very tall peak around 0, suggesting that it may be more susceptible to
being biased by outliers. The spread of changes in emissions is less skewed,
and would make for a better treatment variable. Figure 2 compares the stock
prices (initialized at 1 USD at the start of 2018) between firms that increase
and decrease their Scope 1 and Scope 2 emissions. From first glance, it
appears as though prices are independent of whether firms increase or decrease
their Scope 1 emissions, while firms that increase Scope 2 emissions fare
better against firms that reduce Scope 2 emissions. However, there could be
correlations with several other variables and a more robust empirical strategy
which controls for other variables and firm-specific effects would be required
for an adequate estimation of carbon premia.

(a) Classified by Scope 1 emissions changes

(b) Classified by Scope 2 emissions changes

Figure 2: Stock prices of firms classified by (a) Scope 1 emissions changes,
and (b) Scope 2 changes, for the year 2018-2019.

| Means (standard deviation)
---|---
| 2018 |  | 2019 |  | $\Delta$
| Emission Reducers | Emission Increasers |  | Emission Reducers | Emission Increasers |  | Emission Reducers | Emission Increasers
Log Total Emissions (MT CO2e) | 13.7206 | 13.8380 |  | 13.5977 | 14.0185 |  | -0.1229 | 0.1806
| (2.4297) | (2.3299) |  | (2.4228) | (2.3059) |  | (0.1977) | (0.4387)
Total Liabilities (Billion USD) | 47.4739 | 33.5895 |  | 50.7103 | 34.6071 |  | 3.2364 | 1.0176
| (91.6004) | (52.9324) |  | (97.8752) | (52.4365) |  | (8.9561) | (4.6629)
Total Assets (Billion USD) | 64.9555 | 49.2088 |  | 68.8202 | 50.4200 |  | 3.8646 | 1.2112
| (110.9986) | (75.0238) |  | (117.7552) | (74.9736) |  | (11.5056) | (9.3826)
Free Cash Flow (Billion USD) | 2.2361 | 3.7334 |  | 2.7821 | 3.5361 |  | 0.5460 | -0.1973
| (4.1155) | (8.8737) |  | (7.1775) | (8.7526) |  | (5.5011) | (2.0240)
Revenue (Billion USD) | 32.3890 | 36.8028 |  | 32.7859 | 40.9709 |  | 0.3969 | 4.1680
| (47.6884) | (64.2901) |  | (49.2913) | (68.7920) |  | (7.7173) | (17.8683)
Operating Expenses (Billion USD) | 5.9696 | 5.9229 |  | 6.3440 | 6.4711 |  | 0.3744 | 0.5482
| (10.6697) | (9.3548) |  | (11.2890) | (10.1253) |  | (1.6191) | (1.1386)
EBITDA (Billion USD) | 5.8406 | 7.3910 |  | 5.6491 | 7.7799 |  | -0.1915 | 0.3889
| (8.0379) | (15.0582) |  | (7.3249) | (15.0398) |  | (1.4937) | (2.5904)
Earnings per Share (USD) | 5.8920 | 5.4490 |  | 6.1197 | 5.6683 |  | 0.2277 | 0.3860
| (4.5913) | (4.3460) |  | (5.2450) | (4.3801) |  | (1.2170) | (1.4849)
Volatility | 0.0164 | 0.0174 |  | 0.0150 | 0.0178 |  | -0.0015 | 0.0004
| (0.0044) | (0.0049) |  | (0.0050) | (0.0062) |  | (0.0029) | (0.0038)

Table 2: Financial characteristics of firms classified as Emission Reducers
and Emission Increasers, in 2018 and 2019. Emission Reducers (106 firms) are
defined as all those firms whose total (Scope 1 + Scope 2) emissions do not
increase between 2018 and 2019, while the rest are defined as Emission
Increases (35 firms). Includes 141 firms from 11 sectors. Averages are
unweighted.

Table 2 presents the means and standard deviations of financial
characteristics for the firms in the sample disaggregated by whether a firm’s
emissions increased between 2018 and 2019 (“Emission Increaser”) or whether
its emissions decreased between 2018 and 2019 (“Emission Reducer”). The sample
is split in almost exactly a three-to-one ratio in favor of emission reducers.
This indicates that not only is the market getting cleaner overall, but the
phenomenon is not due to a minority of firms cutting down on emissions by
large amounts, but rather the dominant firm behavior. There are considerable
disparities in the characteristics between Emission Reducers and Emission
Increasers. As of 2018, Emission Reducers have more liabilities, less free
cash flow, less revenue despite marginally higher operating expenses, and as a
consequence lower EBITDA. However, these companies also have more assets,
higher earnings per share, and slightly lower volatility.444We define
volatility to be the standard deviation of the daily log returns over one
fiscal year. Emission reducers have more assets and liabilities, but emission
increases have higher revenue, profits, and cash flow. This suggests that
firms that have high revenue and profits, may be less incentivized to cut down
on their emissions due to high financial performance, although sectoral
differences in Reducers and Increasers could also play an important role.5554
of 11 Energy companies increased emissions, while the numbers were 4(14) for
Materials, 4(20) for Industrials, 3(9) for Consumer Discretionary, 2(7) for
Consumer Staples, 5(20) for Health Care, 0(12) for Financials, 8(22) for
Information Technology, 0(1) for Communication Services, 3(20) for Utilities,
and 2(5) for Real Estate. If existing firm characteristics influenced whether
a firm increased their emissions or not, then a selection bias in assignment
to treatment must be corrected for to study whether there exists a carbon
premium, or not.

##  3 Identification Strategy

###  3.1 Difference-in-differences design

In order to study the effect of GHG emissions on stock prices, we begin with a
two-way fixed effects model (Roth et al., 2023), with the logarithm of returns
$r_{ijt}$ of firm $i$ of sector $j$ at time $t$, as the outcome variable:

$r_{ijt}=\mu_{i}+\alpha_{t}+\theta_{j}+\beta\cdot E_{ijt}+\gamma\cdot
X_{ijt}+\epsilon_{ijt},~{}~{}~{}t=2018,2019$ (1)

Here, $\mu_{i}$ are firm-specific effects, $\alpha_{t}$ are time-specific
effects, $\theta_{j}$ are sector-specific effects, $E_{ijt}$ is the vector of
firms’ Scope 1 and Scope 2 emissions, while $\beta$ is the vector of
coefficients that measures the impact of the firms’ Scope 1 and Scope 2
emissions, and $\gamma_{ijt}$ is the vector of coefficients of the control
variables in vector $X_{ijt}$, and finally, $\epsilon_{ijt}$ is the
idiosyncratic error term. Here, if $\beta$ is estimated as a statistically
significant positive constant, it would imply that a positive carbon premium
exists in the market. On the contrary, if $\beta$ is found to be negative, it
means an equity greenium exists. Note that while we have data on
$r_{ijt},E_{ijt}$, and $X_{ijt}$, we do not have data on $\alpha_{t},\mu_{i},$
and $\theta_{j}$. Since $\mu_{i}$ and $\theta_{j}$, however, are time-
invariant effects, we can eliminate these terms by taking the first difference
of equation 1, and arrive at

$\Delta r_{ijt}=r_{ijt}-r_{ijt-1}=\Delta\alpha_{t}+\beta\cdot\Delta
E_{ijt}+\gamma\cdot\Delta X_{ijt}+\Delta\epsilon_{ijt},~{}~{}~{}t=2019$ (2)

Since we have data for all of $\Delta r_{ijt},\Delta E_{ijt},$ and $\Delta
X_{ijt}$, we can estimate $\beta$ and $\gamma$ using ordinary least squares
(OLS), with $\Delta\alpha_{t}$ as the intercept.

###  3.2 Inverse propensity score weighting

A simple difference-in-differences specification assumes that the assignment
of treatment is “as good as random”, i.e., whether a firm increases or reduces
its emissions, and by how much, can be assumed to be determined randomly.
However, this need not be true. Without additional tests, we cannot ascertain
whether emissions cause an increase in stock prices or whether it is those
firms that bring high returns to investors that are more lax about
environmental governance. While the two are not mutually exclusive, it is
imperative to correct for the latter while estimating the first, i.e., we need
to ensure that there is no selection bias in the assignment of firms to
treatment.

One approach to control for selection bias in assignment to treatment is to
weight observations by the inverse of their propensity scores to receiving the
treatment that they received. The rationale behind controlling for selection
bias assumes that firms are not “treated” at random, i.e., whether a firm
increases or decreases emissions is not random and may depend on the firm’s
existing covariates. Propensity scores are estimated to indicate the
probability of a firm receiving the treatment that it got based on observable
covariates, dividing by which results in a quasi-experimental set-up where
selection bias is eliminated. The intuitive basis for inverse propensity score
weighting is as follows: In a hypothetical case with a binary treatment,
suppose, firms are classified into those that have profited in the preceding
year, and those that have made losses. Suppose a profiting firm is twice as
likely to increase emissions as a loss-making firm, i.e, say, the probability
of inclusion to the treatment group for profit-making firms is 0.5, while that
of loss-making firms is 0.25. Thus, there would be twice as many profit-making
firms that increased emissions, when compared to loss making firms, which
makes the treatment non-random. Thus, by weighting each firm by the inverse of
the propensity score, here, assigning a weight of 2 for profit-making firms,
and a weight of 4 for loss-making firms, both type of firms are equally likely
to be represented in the treatment group. Similarly, in the control group,
profit-making firms would have a higher weight, since they are normally less
likely to be in the control group, in this example.

In our case, the treatment variable is the change in the log emissions (see
equation 2). Thus we have a case of estimating generalized propensity scores
(Hirano and Imbens, 2004) with continuous treatment doses. Following Hirano
and Imbens (2004), we assume a normal distribution for treatment-given
covariates. Let $Z_{ijt}$ be the vector of covariates that confound the
assignment to treatment $\Delta E_{ijt}$, then the coefficients of confounders
in determining the treatment can be estimated from the following equation:

$\Delta E_{ijt}=\beta_{0}^{\prime}+\beta_{1}^{\prime}\cdot
Z_{ijt}+\varepsilon_{ijt}$ (3)

where $\varepsilon$ has the distribution $\mathcal{N}(0,\sigma^{2})$. Having
estimated the coefficients of $\beta_{0}^{\prime}$ and $\beta_{1}^{\prime}$,
we now have an estimate for the expected treatment dose $\hat{\Delta E_{ijt}}$
given by

$\hat{\Delta E_{ijt}}=\beta_{0}^{\prime}+\beta_{1}^{\prime}\cdot Z_{ijt}$ (4)

Using ${\Delta E_{ijt}}$ and $\hat{\Delta E_{ijt}}$, we estimate the
generalized propensity score $p_{i}$ for unit $i$ as

$\hat{p_{i}}=\frac{1}{\sqrt{2\pi\hat{\sigma}^{2}}}\left(-\frac{1}{2\hat{\sigma}^{2}}\left(\Delta
E_{ijt}-\hat{\Delta E_{ijt}}\right)^{2}\right)$ (5)

We can then weight each firm by

$w_{i}=\frac{1}{\hat{p_{i}}}$ (6)

Note that while the propensity-score-weighted-DID regression controls for
selection bias due to observables, we still need to make the unconfoundedness
assumption, i.e., we assume that there are no unobservable characteristics
that induce bias in the treatment that each firm receives.

##  4 Results

###  4.1 What factors cause firms to increase or decrease emissions?

For a simple difference-in-differences analysis, it is important that the
treatment be random, i.e., changes in firms’ carbon emissions would need to be
exogenously determined by firm-level decisions alone independent of other
covariates such as size, profitability, etc., which is unlikely to be true.
Thus, we begin by studying how the changes in emissions of firms depend on
firm characteristics in the pre-treatment period. Table 3 presents the results
of the effects of various 2018 levels of firm characteristics on subsequent
changes in GHG emissions between 2018 and 2019 for both Scope 1 and Scope 2
emissions. We can observe that the change in Scope 2 emissions is not affected
statistically significantly by any covariate. This implies that there may not
be fixed reasons why firms change their Scope 2 emissions, and they may simply
be determined by firm decisions or other unobservable characteristics, which
we cannot test.

| Coefficients (standard errors)
---|---
|  $\Delta$ Log Scope 1 Emissions |  |  $\Delta$ Log Scope 2 Emissions
Intercept | -1.4359* |  | -1.9809
| (0.786) |  | (2.825)
Log Total Liabilities | -0.2635** |  | -0.0762
| (0.102) |  | (0.367)
Log Total Assets | 0.3344*** |  | 0.0936
| (0.116) |  | (0.418)
Free Cash Flow† | -0.0871 |  | 0.3957
| (0.214) |  | (0.769)
Log Revenue | -0.0088 |  | 0.1002
| (0.062) |  | (0.222)
Operating Expenses† | -0.1933 |  | -0.5437
| (0.226) |  | (0.812)
EBITDA† | 0.3299 |  | 0.2967
| (0.242) |  | (0.870)
Earnings per Share $\times 10^{-7}$ | -0.0015 |  | 0.0941
| (0.084) |  | (0.303)
Volatility | 12.4585* |  | 0.6886
| (6.607) |  | (23.742)
Log Scope 1 Emissions | -0.0467*** |  | -0.0600
| (0.014) |  | (0.052)
Log Scope 2 Emissions | 0.0203 |  | -0.0538
| (0.017) |  | (0.062)
Momentum | 0.0819 |  | 0.0227
| (0.141) |  | (0.508)
Sector Fixed Effects | Yes |  | Yes

Table 3: The effect of 2018-levels of firm characteristics on the change in
GHG emissions of firms between 2018 and 2019, over a sample of 141 firms
indexed in S&P500. Heteroskedasticity consistent standard errors. †Variable
scaled by revenue.
***Significant at the 1% level **Significant at the 5% level *Significant at
the 10% level

On the other hand, changes in Scope 1 emissions are strongly correlated with
several pre-treatment characteristics. A firm is likely to observe a more
negative change in Scope 1 emissions if it has fewer assets, more liabilities,
less volatility, and high existing Scope 1 emissions, with all effects being
statistically significant (at least at the 5% level). The last is of
particular relevance because it implies that big polluters are self-correcting
their emissions, although the effect itself is rather small. The EBITDA also
has a positive effect on the change in Scope 1 emissions but is significant
only at the 10% level. However, since we study a relatively small sample of
141 firms, this effect cannot be dismissed entirely, as standard errors can be
expected to be larger in smaller samples. Note that while Scope 2 emissions
reduce by nearly 2% on average (having controlled for observables), which is
almost one and a half times the change for Scope 1 emissions, this effect is
statistically insignificant.

While these results are valuable in isolation as well, they are also important
for the estimation of the effects of GHG emissions on stock returns. Since
changes in Scope 1 emissions are not random, the causal effect cannot merely
be studied through levels regressions, fixed effects regressions, or even
simple difference-in-differences regressions, and the selection bias in
assignment to treatment will need to be controlled for, before effects can be
appropriately estimated.

###  4.2 Is there a carbon premium?

Table 4 presents the results of the unweighted (column 1) and weighted DID
regressions, for the impact of emissions on stock returns, along with some
important control variables, with column 2 presenting the results with Scope 1
emissions as the treatment and column 3 presenting the results with Scope 2
emissions as the treatment. The coefficients presented in column 1 and column
2, respectively, are our preferred estimate of the risk-premium associated
with Scope 1 emissions and Scope 2 emissions. We find that there is a
statistically significant (at the 1% level) risk-premium associated with Scope
1 emissions. For every 1% increase in Scope 1 emissions, there is a 0.076%
surplus return in the shares of a stock. The effect is of a similar size in
the unweighted regression but is only significant at the 10% level. This
signifies that when firms that are less likely to increase emissions actually
increase their emissions, there is a bigger rise in the stock price. This
captures the effect of an unanticipated increase in a firm’s polluting
behavior, leading to a greater pricing of carbon risk in the market for said
firm.

| Coefficients (standard errors)
---|---
| Unweighted DID | p-weighted DID
|  | Weighting: Scope 1 Emissions | Weighting: Scope 2 Emissions
| (1) | (2) | (3)
Intercept | 0.2176*** | 0.2191*** | 0.2453*
| (0.0476) | (0.0466) | (0.1453)
$\Delta$ Log Scope 1 Emissions | 0.0629* | 0.0761*** | 0.1113
| (0.0344) | (0.0245) | (0.1004)
$\Delta$ Log Scope 2 Emissions | -0.0099 | -0.0148 | 0.0114
| (0.0356) | (0.0338) | (0.0110)
$\Delta$ Log Total Liabilities | -0.2230 | -0.3530 | -0.5456
| (0.2334) | (0.2465) | (0.4463)
$\Delta$ Log Total Assets | 0.4574 | 0.6128* | 0.9750
| (0.3144) | (0.3191) | (0.6791)
$\Delta$ Free Cash Flow | -0.1709 | -0.1930 | -0.1378
| (0.1207) | (0.1203) | (0.1657)
$\Delta$ Log Revenue | -0.2775** | -0.2727** | -0.6478
| (0.1175) | (0.1063) | (0.3953)
$\Delta$ Operating Expenses | 0.0318 | 0.1701 | -2.1917
| 0.3961 | (0.3971) | (2.1872)
$\Delta$ EBITDA | 0.3106 | 0.3432 | 1.2074
| (0.3658) | (0.3771) | (1.0857)
$\Delta$ Earnings per Share | 0.0000 | 0.0000 | 0.0000
| (0.0109) | (0.0068) | (0.0200)
$\Delta$ Volatility | 1.0693 | 0.9510 | -9.7905
| (5.6392) | (5.1904) | (16.2962)
Sector Fixed Effects | Yes | Yes | Yes

Table 4: The results for the impact of GHG emissions (carbon premium) and
select control variables on stock returns of firms between 2018 and 2019, over
a sample of 141 firms indexed in S&P500. Heteroskedasticity consistent
standard errors.
***Significant at the 1% level **Significant at the 5% level *Significant at
the 10% level

However, we find that Scope 2 emissions do not produce a statistically
significant risk premium, although the size of the premium is bigger. For
every 1% increase in Scope 2 emissions, there is a 0.111% abnormal excess
return on the stock that cannot be explained by other observable
characteristics. However, due to high variation in the behavior of stocks,
this effect is not significant. It is noteworthy that in the absence of
weighting by propensity scores it may seem that there is a statistically
insignificant equity greenium for Scope 2 stocks which would have led to a
fall in stock prices for a rise in emissions, however, controlling for the
assignment to treatment, the effect is found to be positive.

###  4.3 Is the carbon premium symmetric?

It is possible that the carbon premium exists only for an increase in carbon
emissions, and reductions aren’t thought of to have significantly lead to a
reduction in the price of the risk. To check this, we consider two treatment
variables each for Scope 1 and Scope 2 emissions, allowing for separate carbon
premia for emissions increasers and reducers.

| Coefficients (standard errors)
---|---
| Unweighted DID | p-weighted DID
|  | Weighting: Scope 1 Emissions | Weighting: Scope 2 Emissions
| (1) | (2) | (3)
Intercept | 0.2416*** | 0.2475*** | 0.3150***
| (0.0522) | (0.0540) | (0.1084)
$\Delta$ Log Scope 1 Emissions – Increasers | 0.0657** | 0.0688*** | 0.0897
| (0.0301) | (0.0262) | (0.0717)
$\Delta$ Log Scope 1 Emissions – Decreasers | 0.0175 | 0.1123 | -0.1529
| (0.1426) | (0.1670) | (0.3258)
$\Delta$ Log Scope 2 Emissions – Increasers | -0.0991 | -0.0909 | -0.0806
| (0.1395) | (0.1557) | (0.2626)
$\Delta$ Log Scope 2 Emissions – Decreasers | 0.0025 | 0.0007 | 0.0160
| (0.0216) | (0.0244) | (0.0108)
$\Delta$ Log Total Liabilities | -0.1993 | -0.2951 | -0.2610
| (0.2349) | (0.2436) | (0.3018)
$\Delta$ Log Total Assets | 0.4042 | 0.5069 | 0.4654
| (0.3174) | (0.3273) | (0.4049)
$\Delta$ Free Cash Flow | -0.1785 | -0.2023 | -0.2140
| (0.1200) | (0.1262) | (0.1803)
$\Delta$ Log Revenue | -0.2597** | -0.2406** | -0.5642
| (0.1202) | (0.1100) | (0.3539)
$\Delta$ Operating Expenses | -0.0018 | 0.2736 | -1.5865
| (0.4414) | (0.4356) | (2.0104)
$\Delta$ EBITDA | 0.2771 | 0.3999 | 1.2132
| (0.3863) | (0.3910) | (1.1494)
$\Delta$ Earnings per Share | 0.0000 | 0.0000 | 0.0000
| (0.0074) | (0.0023) | (0.0113)
$\Delta$ Volatility | 2.7817 | 2.1787 | 4.7551
| (5.5513) | (4.9504) | (7.6917)
Sector Fixed Effects | Yes | Yes | Yes

Table 5: The results for the impact of GHG emissions (carbon premium) and
select control variables on stock returns of firms between 2018 and 2019, over
a sample of 141 firms indexed in S&P500. Heteroskedasticity consistent
standard errors.
***Significant at the 1% level **Significant at the 5% level *Significant at
the 10% level

Table 5 presents the results for the estimation of carbon premia with
different estimates for firms that increase and decrease their emissions.
Column 2, our preferred estimate for the carbon premia associated with Scope 1
emissions, shows that carbon risk associated with Scope 1 emissions are only
priced when emissions increase, and that there is no alleviation of carbon
risk when emissions reduce at the rate at which they have, implying that
emissions would have to reduce at a much faster rate in order to reduce risk
for investors. Ideally, we would’ve expected to find a similar effect for
Scope 2 emissions (Column 3), however, no effect on the premium is found on
either increases or decreases in Scope 2 emissions. However, it is difficult
to draw inferences from the estimates of Scope 2 emissions of emissions
increasers as there are only 20 firms of this nature in the sample.

##  5 Conclusion

While several earlier studies have explored the presence of a carbon premium
or equity greenium in financial markets, the causal aspect of the link has
received little attention. In this study, we investigate whether there exists
a causal link between emissions and the stock prices of firms. Using a panel
of 141 stocks from the S&P500 companies, we study if Scope 1 and Scope 2
emissions significantly impact stock prices. To deal with the problems of
reverse causality and selection bias, we use a propensity-score-weighted
difference-in-differences analysis to isolate the impact of emissions on
prices. Our first stage regression finds that several firm characteristics
such as a firms assets, liabilities, volatility and Scope 1 emissions have
statistically significant effects on changes in firms’ emissions. We then use
these estimates to control for selection bias by propensity score weighting in
the seconds stage. By scaling by propensity scores, we effectively balance the
sample such that each firm is equally likely to receive a treatment (changes
in emissions). Using a difference-in-differences regression which controls for
sector, time, and firm-fixed effects, we find that there is a statistically
significant, positive carbon premium associated with Scope 1 emissions.
However, we do not find such an impact for Scope 2 emissions. We hypothesize
that this may be because Scope 1 emissions increase on average, while Scope 2
emissions decrease, and an increase in emissions could lead to an increase in
the price of carbon risk but a reduction at the current rate need not
alleviate the risk. To test this, we estimate the coefficients separately for
firms that increase emissions and those that decrease emissions. We find that
for Scope 1 emissions which has a carbon premium overall, the carbon premium
is largely due to firms that increase emissions and firms that decrease
emissions do not have a carbon premium associated with them, which may confirm
our hypothesis.

Our results agree with the Bolton and Kacperczyk (2021b) in the presence of a
significant positive premium associated with Scope 1 emissions in US markets,
but differ in the context of Scope 2 emissions. The difference may be due to a
more rapid decarbonization of Scope 2 emissions in recent years reflected in
the data used in our analysis compared to the flat or weak increases in Scope
2 emissions in the period considered by Bolton and Kacperczyk (2021b). It is
also possible that our estimate of Scope 2 emissions is not significant due to
a smaller sample of firms. Nevertheless, our results lend support to existing
studies on carbon premia, by eliminating biases typically present in
correlation regressions, by controlling for assignment to treatment.

Our findings reveal key insights about the pricing of carbon risk in financial
markets. First, we find that investors are aware of carbon risk and are
demanding greater short-term returns from polluting firms. Secondly, we find
that investors expect firms to reduce their emissions over time, and thus a
reduction in emissions at the current rate does not alleviate risks associated
with polluting firms. This implies that firms must decarbonize at a much
faster rate to reduce investor risks. On the contrary, increases in firms’s
emissions leads to a substantial increase in the carbon risk which is
demonstrated by the presence of carbon premia for stocks of companies that
increase their emissions.

## References

  * Alessi et al. (2021) Lucia Alessi, Elisa Ossola, and Roberto Panzica.  What greenium matters in the stock market? the role of greenhouse gas emissions and environmental disclosures.  _Journal of Financial Stability_ , 54:100869, 2021.
  * Andersson et al. (2016) Mats Andersson, Patrick Bolton, and Frédéric Samama.  Hedging climate risk.  _Financial Analysts Journal_ , 72(3):13–32, 2016\.
  * Aswani et al. (2022) Jitendra Aswani, Aneesh Raghunandan, and Shivaram Rajgopal.  Are carbon emissions associated with stock returns?  _Columbia Business School Research Paper Forthcoming_ , 2022.
  * Bernardini et al. (2021) Enrico Bernardini, Johnny Di Giampaolo, Ivan Faiella, and Riccardo Poli.  The impact of carbon risk on stock returns: evidence from the european electric utilities.  _Journal of Sustainable Finance & Investment_, 11(1):1–26, 2021.
  * Bolton and Kacperczyk (2021a) Patrick Bolton and Marcin Kacperczyk.  Global pricing of carbon-transition risk.  Technical report, National Bureau of Economic Research, 2021a.
  * Bolton and Kacperczyk (2021b) Patrick Bolton and Marcin Kacperczyk.  Do investors care about carbon risk?  _Journal of financial economics_ , 142(2):517–549, 2021b.
  * Brunborg and Haldorsen (2021) Thor-Bjørn Brunborg and Per Kristian Haldorsen.  The impact of carbon emissions on us stock returns: An empirical study of how carbon intensity impacts financial performance.  Master’s thesis, 2021.
  * Chakrabarty and Nag (2023) Siddhartha P Chakrabarty and Suryadeepto Nag.  Risk measures and portfolio analysis in the paradigm of climate finance: a review.  _SN Business & Economics_, 3(3):69, 2023.
  * Cucchiella et al. (2017) F Cucchiella, M Gastaldi, and M Miliacca.  The management of greenhouse gas emissions and its effects on firm performance.  _Journal of Cleaner production_ , 167:1387–1400, 2017.
  * Curtin et al. (2019) J Curtin, C McInerney, B Ó Gallachóir, C Hickey, P Deane, and P Deeney.  Quantifying stranding risk for fossil fuel assets and implications for renewable energy investment: a review of the literature.  _Renewable and Sustainable Energy Reviews_ , 116:109402, 2019.
  * Del Prado (2019) I Del Prado.  In the view of a carbon tax: How carbon emission and climate change action are priced in us stocks.  2019\.
  * Di Cosmo and Hyland (2013) Valeria Di Cosmo and Marie Hyland.  Carbon tax scenarios and their effects on the irish energy sector.  _Energy Policy_ , 59:404–414, 2013.
  * Fawcett (2004) Tina Fawcett.  Carbon rationing and personal energy use.  _Energy & Environment_, 15(6):1067–1083, 2004\.
  * Görgen et al. (2020) Maximilian Görgen, Andrea Jacob, Martin Nerlinger, Ryan Riordan, Martin Rohleder, and Marco Wilkens.  Carbon risk.  _Available at SSRN 2930897_ , 2020.
  * Griffin (2017) Paul Griffin.  CDP carbon majors report 2017.  Technical report, Carbon Disclosure Project, 2017.
  * Gurvich and Creamer (2022) Alex Gurvich and Germán G Creamer.  Carbon risk factor framework.  _The Journal of Portfolio Management_ , 48(10):148–164, 2022.
  * Hamilton and Cameron (1994) Kirk Hamilton and Grant Cameron.  Simulating the distributional effects of a canadian carbon tax.  _Canadian Public Policy/Analyse de Politiques_ , pages 385–399, 1994\.
  * Hirano and Imbens (2004) Keisuke Hirano and Guido W Imbens.  The propensity score with continuous treatments.  _Applied Bayesian modeling and causal inference from incomplete-data perspectives_ , 226164:73–84, 2004.
  * Howell (2012) Rachel A Howell.  Living with a carbon allowance: The experiences of carbon rationing action groups and implications for policy.  _Energy Policy_ , 41:250–258, 2012.
  * Ji et al. (2021) Xiangfeng Ji, Yusong Zhang, Nawazish Mirza, Muhammad Umar, and Syed Kumail Abbas Rizvi.  The impact of carbon neutrality on the investment performance: Evidence from the equity mutual funds in brics.  _Journal of Environmental Management_ , 297:113228, 2021\.
  * Kim et al. (2015) Yeon-Bok Kim, Hyoung Tae An, and Jong Dae Kim.  The effect of carbon risk on the cost of equity capital.  _Journal of Cleaner Production_ , 93:279–287, 2015.
  * Krueger et al. (2020) Philipp Krueger, Zacharias Sautner, and Laura T Starks.  The importance of climate risks for institutional investors.  _The Review of Financial Studies_ , 33(3):1067–1111, 2020.
  * Kuang and Liang (2020) Huan Kuang and Bing Liang.  Carbon risk exposure in the mutual fund industry.  _Available at SSRN 3750244_ , 2020.
  * Lo et al. (2022) Andrew W Lo, Ruixun Zhang, and Chaoyi Zhao.  Measuring and optimizing the risk and reward of green portfolios.  _The Journal of Impact and ESG Investing_ , 2022.
  * Maddox (2018) Grace Emily Maddox.  _Emissions trading, carbon pricing, and the impact on carbon-producing firms: a study of phase III of the European Union Emissions Trading Scheme_.  PhD thesis, Massey University, 2018.
  * Mazzarano et al. (2021) Matteo Mazzarano, Gianni Guastella, Stefano Pareglio, Anastasios Xepapadeas, et al.  Carbon boards and transition risk: Explicit and implicit exposure implications for total stock returns and dividend payouts.  Technical report, JSTOR, 2021.
  * Meyer (2022) Fabian Alexander Meyer.  Carbon risk in european equity returns.  _Junior Management Science_ , 7(2):429–454, 2022\.
  * Muñoz (2021) Fernando Muñoz.  Carbon-intensive industries in socially responsible mutual funds’ portfolios.  _International Review of Financial Analysis_ , 75:101740, 2021.
  * Nag et al. (2023) Suryadeepto Nag, Siddhartha P Chakrabarty, and Sankarshan Basu.  Single event transition risk: A measure for long term carbon exposure.  _MethodsX_ , 10:102001, 2023.
  * Oestreich and Tsiakas (2015) A Marcel Oestreich and Ilias Tsiakas.  Carbon emissions and stock returns: Evidence from the eu emissions trading scheme.  _Journal of Banking & Finance_, 58:294–308, 2015.
  * Oswald et al. (2023) Yannick Oswald, Joel Millward-Hopkins, Julia K Steinberger, Anne Owen, and Diana Ivanova.  Luxury-focused carbon taxation improves fairness of climate policy.  _One Earth_ , 2023.
  * Park et al. (2022) Dojoon Park, Jiyoon Lee, and Hyejin Park.  The asset-pricing implications of carbon risk in korea.  _Available at SSRN 4257054_ , 2022.
  * Pástor et al. (2021) L’uboš Pástor, Robert F Stambaugh, and Lucian A Taylor.  Sustainable investing in equilibrium.  _Journal of Financial Economics_ , 142(2):550–571, 2021.
  * Pastor et al. (2022) Lubos Pastor, Robert F Stambaugh, and Lucian A Taylor.  Dissecting green returns.  _Journal of Financial Economics_ , 146(2):403–424, 2022.
  * Ravina and Hentati Kaffel (2019) Alessandro Ravina and Rania Hentati Kaffel.  The impact of low-carbon policy on stock returns.  _Available at SSRN 3444168_ , 2019.
  * Roncalli et al. (2020) Théo Roncalli, Théo Le Guenedal, Frédéric Lepetit, Thierry Roncalli, and Takaya Sekine.  Measuring and managing carbon risk in investment portfolios.  _arXiv preprint arXiv:2008.13198_ , 2020.
  * Roncalli et al. (2021) Théo Roncalli, Théo Le Guenedal, Frédéric Lepetit, Thierry Roncalli, and Takaya Sekine.  The market measure of carbon risk and its impact on the minimum variance portfolio.  _The Journal of Portfolio Management_ , 47(9):54–68, 2021.
  * Roth et al. (2023) Jonathan Roth, Pedro HC Sant’Anna, Alyssa Bilinski, and John Poe.  What’s trending in difference-in-differences? a synthesis of the recent econometrics literature.  _Journal of Econometrics_ , 2023.
  * Rutherford (1992) Thomas Rutherford.  The welfare effects of fossil carbon restrictions: results from a recursively dynamic trade model.  1992\.
  * Santi and Moretti (2021) Caterina Santi and Angelo Moretti.  Carbon risk premium and worries about climate change.  _Available at SSRN 3942738_ , 2021.
  * Semieniuk et al. (2021) Gregor Semieniuk, Emanuele Campiglio, Jean-Francois Mercure, Ulrich Volz, and Neil R Edwards.  Low-carbon transition risks for finance.  _Wiley Interdisciplinary Reviews: Climate Change_ , 12(1):e678, 2021.
  * Tucker (2010) Jenny Wu Tucker.  Selection bias and econometric remedies in accounting and finance research.  _Journal of Accounting Literature_ , 29:31–57, 2010.
  * Vymyatnina and Chernykh (2022) Yulia Vymyatnina and Alexander Chernykh.  Green factor influence on the yield of stocks and bonds in the russian financial market.  _Available at SSRN_ , 2022.
  * Walkshäusl (2021) Christian Walkshäusl.  Carbon momentum.  _The Journal of Impact and ESG Investing_ , 2(1):33–46, 2021.
  * Wen et al. (2020) Fenghua Wen, Nan Wu, and Xu Gong.  China’s carbon emissions trading and stock returns.  _Energy Economics_ , 86:104627, 2020.
  * Witkowski et al. (2021) Pawel Witkowski, Adam Adamczyk, and Slawomir Franek.  Does carbon risk matter? evidence of carbon premium in eu energy-intensive companies.  _Energies_ , 14(7):1855, 2021.
  * Zhang et al. (2018) Miao Zhang, Russell B Gregory-Allen, et al.  Carbon emissions and stock returns: Evidence from the chinese pilot emissions trading scheme.  _Theoretical Economics Letters_ , 8(11):2082, 2018\.
  * Zhang (2022) Shaojun Zhang.  Do investors care about carbon risk? a global perspective.  _A Global Perspective (May 17, 2022)_ , 2022.
  * Zhang and Wei (2010) Yue-Jun Zhang and Yi-Ming Wei.  An overview of current research on eu ets: Evidence from its operating mechanism and economic effect.  _Applied Energy_ , 87(6):1804–1814, 2010.


# Prime Factorization Using Quantum Variational Imaginary Time Evolution

Raja Selvarajan  Department of Chemistry, Department of Physics and Astronomy,
and Purdue Quantum Science and Engineering Institute, Purdue University, West
Lafayette, IN 47907 USA  Vivek Dixit  Department of Chemistry, Department of
Physics and Astronomy, and Purdue Quantum Science and Engineering Institute,
Purdue University, West Lafayette, IN 47907 USA  Xingshan Cui  ‘Department of
Mathematics, Department of Physics and Astronomy, and Purdue Quantum Science
and Engineering Institute, Purdue University, West Lafayette, IN 47907 USA’
Quantum Science Center Oak Ridge National Laboratory, Oak Ridge, TN  Travis S.
Humble This manuscript has been authored by UT-Battelle, LLC, under contract
DE-AC05-00OR22725 with the US Department of Energy (DOE). The US government
retains and the publisher, by accepting the article for publication,
acknowledges that the US government retains a nonexclusive, paid-up,
irrevocable, worldwide license to publish or reproduce the published form of
this manuscript, or allow others to do so, for US government purposes. DOE
will provide public access to these results of federally sponsored research in
accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-
public-access-plan). Quantum Science Center Oak Ridge National Laboratory, Oak
Ridge, TN  Sabre Kais  Department of Chemistry, Department of Physics and
Astronomy, and Purdue Quantum Science and Engineering Institute, Purdue
University, West Lafayette, IN 47907 USA  Quantum Science Center Oak Ridge
National Laboratory, Oak Ridge, TN<EMAIL_ADDRESS>
###### Abstract

The road to computing on quantum devices has been accelerated by the promises
that come from using Shor’s algorithm to reduce the complexity of prime
factorization. However, this promise hast not yet been realized due to noisy
qubits and lack of robust error correction schemes. Here we explore a
promising, alternative method for prime factorization that uses well-
established techniques from variational imaginary time evolution. We create a
Hamiltonian whose ground state encodes the solution to the problem and use
variational techniques to evolve a state iteratively towards these prime
factors. We show that the number of circuits evaluated in each iteration
scales as $O(n^{5}d)$, where $n$ is the bit-length of the number to be
factorized and $d$ is the depth of the circuit. We use a single layer of
entangling gates to factorize several numbers represented using 7, 8, and
9-qubit Hamiltonians. We also verify the method’s performance by implementing
it on the IBMQ Lima hardware.

##  1 Introduction

Quantum computation is likely to revolutionize how computation is performed in
the field of science, engineering, and finance. Computations are performed on
quantum states that make use of superposition and entanglement to allow for
speedups. Future potential applications include cryptography [1], search
problems [2], simulation of quantum systems [3], quantum annealing [4],
machine learning [5], computation biology [6], quantum materials [7], and
problems in optimization [8]. Refer [9] for an extensive introduction into the
field of quantum computation. Major efforts toward scaling the current
algorithm focuses on developing error correction and mitigation schemes, as
well as designing operations that make use of fewer ancilla qubits and gate
operations [10]. In line with the current major developments, we investigate a
more practical near-term scheme for prime factorization that is likely to
achieve good results on noisy qubits.

Prime factorization involves expressing a composite number as the product of
its prime factors. For generic large numbers that lack any structure, the
quadratic sieve is the most commonly used technique. This can be
computationally expensive and is exploited in RSA cryptography to guarantee
information security over networks. From Shor’s [11] work on period finding,
it was shown that one could exploit the quantum Fourier transform to compute
factors in steps that scaled polynomial in the number of bits. Subsequently,
Vandersypen et al. [12] realized this experimentally by factorizing $N=15$
using spin-$1/2$ nuclei as qubits and then Lucero et al. [13] by using
superconducting qubits. A simplified version of Shor’s algorithm was worked
out for products of Fermat primes (3, 5, 17, 257 and 65537) by Geller et al.
[14]. Jian et al. [15] proposed an alternative method to compute factors by
solving an optimization problem using quantum annealing demonstrated on the
D-Wave quantum annealer. The largest experimental realization of general
method factorization schemes includes Shor’s algorithm to factorize 21 [13]
and optimization using the D-Wave quantum annealer to factorize 223357 [15].
In addition, large numbers with specific properties have also been factorized
by exploiting structure contained in the number. In Ref. [16], the authors use
a multiplication table to factorize 56153 using only 4 qubits. Despite being
able to factorize large numbers with relatively few qubits by exploiting the
structure of the number, these methods do not reveal the power of quantum
computers against generic numbers. Studies have not been made with respect to
convergence on the solution with increasing the number of qubits.

In this paper we explore how one could use imaginary time evolution to
factorize numbers with relatively higher probability. Imaginary time evolution
has long been used as a theoretical tool in physics to compute ground-state
wave-functions [17]. Shingu et al. [18] show how using imaginary time
evolution one could train a Boltzmann network efficiently, computing the exact
model independent term in the training. Recently, Motta et al. [19] exploits
the Quantum Imaginary Time Evolution (QITE) algorithm to determine eigenstates
and thermal states on a quantum computer. McArdle et al. [20] showed how one
could make use of variational circuits to create states that represent the
dynamics of the imaginary time evolution. They use it to compute the ground
state wavefunction of hydrogen and lithium hydride, while Yeter-Aydeniz et al.
[21] demonstrate the QITE algorithm as a quantum computing benchmark for
computational chemistry methods.

Herein, we develop an optimization function using the method referenced at
[22] and then use it as a Hamiltonian to perform imaginary time evolution on a
uniform superposition of all possible considered solutions to the
factorization problem. We employ variational circuits to prepare a quantum
state that encodes the solution and classically train all parameters.
Simulations using Python Numpy packages and IBM-QASM are used to verify the
performance. The robustness of the techniques against noise is verified on
IBMQx2 hardware for up to 5 qubits allowing factorization of numbers up to 91
with a probability over $73\%$. To the best of our knowledge this is the
largest number that has been factorized on a quantum circuit using a general
purpose algorithm that does not exploit any structural properties of the
number. Each iteration involves evaluating a number of circuits that scales as
$O(n^{5}d)$ where $n$ is the bit-string length of the number to be factorized
and $d$ is the circuit depth. A few reasons make factorization an ideal
candidate to be solved using Imaginary time evolution with variational
circuits. The number of terms in the Hamiltonian expansion scales polynomial
with respect to the number of qubits. The amplitude of coefficients in the
state vector is real, which simplifies the parameterization of the landscape
to be explored and updates to be made. Each iteration only needs to amplify
the amplitude of the solution rather than exactly simulate the imaginary time
evolution. These factors make factorization using the QITE algorithm a good
candidate to demonstrate the practicality of quantum circuits in near term
devices in the NISQ [3] era.

##  2 Quantum Imaginary Time Evolution

The Schrödinger equation for a closed system describes dynamics according to
some Hamiltonian that governs it. Allowing for time to be complex, we are able
to create thermal states of specific temperature starting from a maximally
mixed configuration, i.e, $\rho_{T=1/\tau}=e^{-H\tau}/$Tr$[e^{-H\tau}]$.
Preparing a system at low temperature or alternatively letting the system
evolve to large imaginary time values we can more often sample the ground
state configuration for a given Hamiltonian. We shall make use of this
property to encode the required solution of the problem we intend to solve in
the later sections.

The quantum state we intend to prepare using QITE is given by,

$\ket{\psi(\tau)}=N(\tau)e^{-H\tau}\ket{\psi(0)}$ (1)

where $N(\tau)=1/\sqrt{\bra{\psi(0)}e^{-2H\tau}\ket{\psi(0)}}$ is a
normalization factor. Equation 1 satisfies the Wick rotated Schrodinger
equation

$\frac{\partial\ket{\psi(\tau)}}{\partial\tau}=-(H-E_{\tau})\ket{\psi(\tau)}$
(2)

where $E_{\tau}=\bra{\psi(\tau)}H\ket{\psi(\tau)}$.

Let $\ket{\phi(\tau)}$ be a state that is prepared by applying a series of
unitary gates as follows

$\ket{\phi(\tau)}=U_{N}(\theta_{N})U_{N-1}(\theta_{N-1})U_{N-2}(\theta_{N-2})U_{N-3}(\theta_{N-3})...U_{1}(\theta_{1})\ket{0}$
(3)

We choose the initial parameters to create the state on which the evolution is
to be performed. We demand that equation (3) is an approximate solution to
equation (2) by demanding the norm of the variation to vanish, i.e,

$\delta\left\lVert\left(\frac{\partial}{\partial\tau}+H-E_{\tau}\right)\ket{\phi(\tau)}\right\rVert=0$
(4)

Solving the above equation using McLachlan’s variational principle (details
expanded in the Appendix) we obtain

$\sum A_{ij}\dot{\theta_{j}}=C_{i}$ (5)

where

$\begin{split}A_{ij}=\frac{\partial\bra{\phi(\tau)}}{\partial\theta_{i}}{\frac{\partial\ket{\phi(\tau)}}{\partial\theta_{j}}}\\\
C_{i}=-\frac{\partial\bra{\phi(\tau)}}{\partial\theta_{i}}H\ket{\phi(\tau)}\end{split}$
(6)

The derivatives on the state vectors with respect to the parameters of the
circuit can be simplified into a sum of other basic circuit functions. For
instance, the derivative of the $R_{x}$, $R_{y}$ and $R_{z}$ rotation gates
with respect to their parameter amounts to a mere shift of the parameter by an
angle of $\pi$ radians. We can further express $H$ as a sum of a string of
tensor products of Pauli operators allowing both $A$ and $C$ to be computed as
the sum of terms that take the form $a\Re(e^{i\phi}\bra{0}U\ket{0})$, where
$a$ is some real number and $\phi$ is a real phase. Terms of this form can be
efficiently computed on quantum hardware using the Hadamard test.

Evaluating the gradients using equation (5), we can then proceed to use
gradient descent to update the parameters of the circuit as follows

$\vec{\theta}(\tau+\delta\tau)=\vec{\theta}(\tau)+\dot{\vec{\theta}}(\tau)\delta\tau$
(7)

where $\dot{\vec{\theta}}(\tau)=A^{-1}(\tau)C(\tau)$. Using sufficient layers
to the variational circuit and small $\delta\tau$ time iterations for updating
our parameters, we can prepare states that closely emulate the imaginary time
evolution on the state $\ket{\psi(0)}$. We would like to note that the circuit
prepares the final state from a fixed starting state rather than from the
state at $\delta t$ time step before evolution as this would correspond to a
non-unitary transformation on the Hilbert space of qubits.

##  3 Prime Factorization to Binary Optimization

The problem of prime factorization involves expressing a number as a product
of its primes. Here we shall focus specifically on biprimes, product of $2$
prime numbers. Shor’s algorithm aims to solve this problem by reducing it to a
period finding problem of the function $a^{x}$ mod $N$, where $a$ is any
random number and $N$ is the number to be factorized. This algorithm uses a
quantum subroutine that exploits the power of quantum Fourier transform as a
part of a modified phase estimation algorithm. The largest number factorized
using Shor’s algorithm is 21. Being highly sensitive to discrete errors,
Shor’s algorithm fails to scale without viable error correction schemes [23].

Let the number $N$ be given by a product of $2$ prime numbers $p$ and $q$. We
model the problem of finding $p$ and $q$ given $N$ as an optimization problem
cast as the following cost function

$H(p,q)=(N-p\times q)^{2}$ (8)

Note that in equation 8, $H(p,q)$ has global minima when $p$,$q$ factorize $N$
with a minimum value of 0. The cost function can thus be treated as a
Hamiltonian whose ground state encodes the solution to our problem. We express
both $p$ and $q$ as a binary string summation which will allow us to transfer
to a qubit representation space. Assuming $p$ has an $m$ bit representation
and $q$ has an $l$ bit representation, we get

$\displaystyle p=2^{m-1}p_{m-1}+2^{m-2}p_{m-2}\ldots 2^{0}q_{0}$ (9)
$\displaystyle q=2^{l-1}q_{l-1}+2^{l-2}q_{l-2}\ldots 2^{0}q_{0}$

Note that equation 8 has a trivial solution $N=N\times 1$. To exclude this
solution, we can restrict the search space of our solution by using a lower
bound. The smallest prime factor we are looking for is greater than $2$ and
less than $N/2$. Thus using $N-1$ qubits to represent $p$ and $q$ we can
discard the trivial solution from showing up in our simulation. Since both $p$
and $q$ are prime numbers, we are allowed to set the bit $0$ to be 1 in either
expression making it indivisible by 2. To move to a scaled spin representation
we transform all the binary variables $b_{i}$ that represents either the
$p_{i}$’s or $q_{i}$’s to $b_{i}=(s_{i}+1)/2$. In this representation, $s_{i}$
takes the value $\pm 1$ such that $s_{i}=1$ maps to $b_{i}=1$ and $s_{i}=-1$
maps to $b_{i}=0$. We thus convert a factorization problem into an
optimization problem over the spin variables.

Minimizing this optimization function shall result in the solution we seek. We
use imaginary time evolution to evolve a uniform superposition over all
possible solutions. The circuit parameters are updated according to equation
(5) and (7). With every iteration that is performed, the output state evolves
through a time step of $\delta\tau$. After sufficiently large enough time $T$,
we expect the ground state with the least cost to prepare. For $\left\lVert
H\right\rVert_{1}T>>1$ we get,

$e^{-HT}\ket{++..+}\approx\ket{\tilde{p}}\ket{\tilde{q}}$ (10)

where $\tilde{p}$ and $\tilde{q}$ represent the bit string solution of the
binary representation of the numbers that multiply to give $N$. The evolution
is achieved by using a variational circuit ansatz to prepares states that
mimic the evolution of the system from a given starting state (as described in
the previous section). We show in the Appendix that the computational cost in
performing one iteration of QITE is $O(n^{5}d)$, where $n$ is the bit-length
of the number to be factored and $d$ is the depth of the circuit.

##  4 Demonstrating an Example

Here we show how to factorize $N=15$ as

$15=5\times 3$ (11)

Expressing the factors 5 and 3 as binary variables and setting $p_{0}=1$,
$q_{0}=1$ yields

$\displaystyle 5=4x_{1}+2x_{0}+1$ (12) $\displaystyle 3=2x_{2}+1$

with the correct results defined by $x_{0}=0,x_{1}=1,x_{2}=1$. The
corresponding cost function $H(p,q)$ is cast as

$H(x_{0},x_{1},x_{2})=[15-(4x_{1}+2x_{0}+1)(2x_{2}+1)]^{2}$ (13)

Expanding the cost function and making use of $x_{i}^{2}=x_{i}$ yields

$H(x_{0},x_{1},x_{2})=196-52x_{2}-52x_{0}-56x_{2}x_{0}-96x_{1}-48x_{2}x_{1}+16x_{0}x_{1}+128x_{0}x_{1}x_{2}$
(14)

Mapping the binary variables to spin variables yields

$H(s_{0},s_{1},s_{2})=90-36s_{2}-40s_{1}-20s_{0}+2s_{2}s_{0}+4s_{2}s_{0}+20s_{0}s_{1}+16s_{0}s_{1}s_{2}$
(15)

The classical binary variables $s_{i}$ in the Hamiltonian is replaced by the
Pauli-z spin operator to get a quantum Hamiltonian Operator that is then acted
upon the full superposition state. Using a time step $\delta\tau=0.1$ , we
make $10$ iterations of QITE to reach $T=1$. With a probability of greater
than $90\%$, we obtain upon measurement

$e^{-HT}\ket{+++}\approx\ket{110}$ (16)

where the above expression is written in the computational basis. The output
equation (16) when mapped back to the binary variables gives

$x_{0}=0,x_{1}=1,x_{2}=1$ (17)

We thus obtain the factors of 3 and 5 as expected. In our simulations, we
shall make use of no more than the number of qubits needed to represent our
solution in the optimization. This is to limit the numbers of qubits used,
maximize the computational efficiency and reduce the accumulated errors.

##  5 Reformulating using QFI and Gradients

We generate the output quantum state prepared using QITE with the ansatz shown
in figure 1. Notice that the circuit consists of $R_{Y}$ rotation gates
applied to each qubit followed by a layer of $CNOT$ gates that help with the
entanglement of the qubits. Using only $R_{Y}$ gates ensures that we maintain
real amplitudes for the state when expressed in the computation basis. This
helps with closely following the dynamics of QITE without introducing any
nontrivial phase values.

Figure 1: Variational circuit used to prepare the quantum state for imaginary
time evolution. The parameters have been randomly initialized

We show how for such an ansatz, equation 5 can be re-expressed in terms of the
quantum Fisher information (QFI) and the gradient of Hamiltonian expectation
with respect to the current state. In classical probability, the Fisher
information characterizes how much a probability distribution varies by
changing a parameter that characterizes a distribution. This can be
generalized and extended to talk about quantum states where one characterizes
how much a state changes with respect to the parameter that governs it. Refer
to [24] for a brief introduction to quantum Fisher information. Given a
quantum state $\ket{\psi(\theta_{1},\theta_{2}...\theta_{n})}$, where
$\theta_{i}$ represents the governing parameters, the QFI is given by

$F_{ij}=4\Re(\braket{\psi_{i}}{\psi_{j}}-\braket{\psi_{i}}{\psi}\braket{\psi}{\psi_{j}})$
(18)

where $\ket{\psi_{j}}=\frac{d\ket{\psi(\theta_{j})}}{d\theta}$ . For the
ansatz being used to generate the state in QITE, we note that
$\braket{\psi_{i}}{\psi}=\Re(\braket{\psi_{i}}{\psi})$. This results in the
second term vanishing due to the constant normalization of the state i.e,
$\braket{\psi}{\psi}=1$. Hence we get $F_{ij}=4M_{ij}$ where
$M_{ij}=\Re(\braket{\psi_{i}(\theta(\tau))}{\psi_{j}(\theta(\tau))})$, cf. Eq.
(26) in Appendix A.

Furthermore, since the Hamiltonian is an Hermitian operator, we can express
the right hand side of equation (5) as follows

$C_{i}=\Re(\bra{\psi}H\ket{\psi_{i}})=\frac{d}{d\theta_{i}}(\bra{\psi}H\ket{\psi})/2$
(19)

We thus obtain
$\sum\limits_{j}F_{ij}\dot{\theta_{j}}=-2\frac{d}{d\theta_{i}}\braket{H}_{\psi}$.
Both the quantum Fisher information of a given circuit and the gradient of an
observable with respect to a state parameter can be directly computed in the
IBM Qiskit Aqua framework.

##  6 Results on Simulation

(a) Amplitude of correct solution for various 5-qubit examples.

(b) Amplitude of correct solution for various 7-qubit examples.

(c) Amplitude of correct solution for various 8-qubit examples.

(d) Amplitude of correct solution for various 9-qubit examples.

Figure 2: Numerical simulations of QITE for 7-, 8-, and 9-qubit factorization
examples. The vertical axis indicates the amplitude of the solution in the
computational basis and the horizontal axis indicates the number of iterations
made. The curves in each subplot show how the amplitude of the solution
corresponding to various numbers increases with iterations. The amplitude
threshold has been set to 0.85.

Smaller numbers were factorized on Qiskit using the QASM simulator. For 8
qubits and more, the simulations were run using the Numpy package. The
amplitude threshold for the correct solution has been set to $0.85$, this
amounts to a probability of $73\%$ to get the right solution when measured in
the computational basis. We have restricted our ansatz to a single layer in
addition to the base layer, in light of exploring shallow circuits for the
near term quantum computers. The initial circuit parameters have been randomly
sampled to help avoid valleys during the training and allow for faster
convergence with shallow circuits employed. We have noticed that it also helps
in suppressing one of the solutions in the case of symmetrically equivalent
solutions (for instance, 391 can be factorized as $17\times 23$ and $23\times
17$), allowing one solution to quickly reach the threshold.Figure 2 plots the
amplitude of the solution in the computational basis against the number of
iterations for several values of $N$ that are represented using 5, 7, 8, and
9-qubit Hamiltonians. Note the amplitude of the solution to the factorization
keeps raising with increasing the number of iterations. The figure also shows
that the average number of iterations required to factorize increases with the
number of qubits though we have not been able to analytically bound it. A few
instances involve the amplitude of the solution flattening before continuing
to raise any further, making it significantly slow. In table 1, we note that
despite some numbers being close to each other (for instance 1007 and 1081),
the number of iterations for convergence to the solution seem to be widely
varying. This can be attributed to the differences in bit string
representation and thus the created cost function, and also to the fact that
the parameters are being randomly initialized. Results on actual hardware for
smaller instances have been presented in section 8.

Number factorized | Number of iterations
---|---
55 | 31
65 | 23
77 | 21
91 | 23
187 | 30
209 | 48
221 | 61
247 | 77
253 | 72
299 | 58
319 | 30
341 | 57
377 | 84
403 | 36
323 | 81
391 | 89
437 | 48
493 | 80

Number factorized | Number of iterations
---|---
527 | 47
551 | 77
589 | 38
629 | 132
697 | 59
703 | 31
713 | 39
731 | 56
799 | 86
899 | 96
1007 | 21
1037 | 155
1081 | 112
1159 | 103
1247 | 112
1457 | 128
1643 | 88
1829 | 130

Table 1:  Table lists the numbers that were factorized alongside the number of
iterations for convergence

##  7 Results on IBM Hardware

We have made use of the publicly available IBM hardware to simulate
factorization for up to 5 qubits. We start the output state with a uniform
superposition over all possible solutions. This is straightforwardly
instantiated by setting the parameters of $R_{y}$ gates in the base layer to
be $\frac{\pi}{2}$ and the remaining parameters to be all $0$. We have avoided
invoking the Hadamard test to compute the $M$ matrix and $C$ vector as this
resulted in significant error accumulation. Instead, we have reformulated the
computation as indicated in section 5 and have made use of built-in Qiskit
modules. The Gradient module with parameter shift method has been used to
compute the $C$ vector, while the QFI module with overlapping block diagonal
method has been used to compute the $M$ matrix. Please refer to the Qiskit
source documentation [25] for details on the implementation of these modules.
Figure 3 shows the results of factorizing the numbers $N=$ 55, 65, 77 and 91
on IBMQ-lima hardware that supports 5 qubits. Unlike the simulations that
hardly show any oscillation in the amplitude, we see oscillations being
introduced due to noise when run on the hardware. The convergence of the
solution in the presence of hardware noise makes this method a suitable
candidate for solving factorization and other similarly framed optimization
problems in the near term quantum computers.

Figure 3: 5-qubit factorization on IBMQ-lima. The vertical axis indicates the
amplitude of the solution in the computational basis and the horizontal axis
indicates the number of iterations made. The curves in each subplot, show how
the amplitude of the solution corresponding to various numbers increase with
iterations. The amplitude threshold has been set to 0.85

##  8 Conclusions

We have shown how imaginary time evolution can be used to perform optimization
and we have demonstrated this for the case of prime factorization. We have
shown numerical adn experimental results from the factorization of several
numbers represented by 7, 8, and 9-qubit Hamiltonians. We have shown that
imaginary time evolution significantly populates the probability of measuring
the correct solution to the factorization problem when run sufficiently long.
The observed performance of the method on the IBM-lima hardware reveals a
robustness towards noise and works as a proof of concept in the case of a
limited number of qubits. Methods on how to bound the number of iterations
analytically can broaden the scope for this method on near-term quantum
computers for very large numbers.

##  9 Acknowledgement

We would like to thank Dr. Mario Motta and Dr. Stefan Woerner of IBM for
useful discussion. We would like to acknowledge the support from National
Science Foundation under award number 1955907. This material is based upon
work supported by the U.S. Department of Energy, Office of Science, National
Quantum Information Science Research Centers, Quantum Science Center.

##  10 Author contributions statement

S.K., T.H. and X.C. designed the research problem, R.S. performed hardware
testing and wrote the paper, R.S and V.D coded for simulations using Python
Numpy. All authors helped in revising the paper draft.

## References

  * [1] Daniel J. Bernstein.  Introduction to post-quantum cryptography, pages 1–14.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2009.
  * [2] Andris Ambainis.  Quantum search algorithms.  ACM SIGACT News, 35:22–35, 05 2005.
  * [3] John Preskill.  Quantum computing in the NISQ era and beyond.  Quantum, 2:79, Aug 2018.
  * [4] Vivek Dixit, Raja Selvarajan, Muhammad A. Alam, Travis S. Humble, and Sabre Kais.  Training and classification using a restricted Boltzmann machine on the D-Wave 2000Q.  preprint arxiv:2005.03247, 2020.
  * [5] Vivek Dixit, Raja Selvarajan, Tamer Aldwairi, Yaroslav Koshka, Mark A. Novotny, Travis S. Humble, Muhammad A. Alam, and Sabre Kais.  Training a quantum annealing based restricted boltzmann machine on cybersecurity data.  IEEE Transactions on Emerging Topics in Computational Intelligence, pages 1–12, 2021.
  * [6] Carlos Outeiral, Martin Strahm, Jiye Shi, Garrett M. Morris, Simon C. Benjamin, and Charlotte M. Deane.  The prospects of quantum computing in computational molecular biology.  WIREs Computational Molecular Science, 11(1):e1481, 2021.
  * [7] Shree Hari Sureshbabu, Manas Sajjan, Sangchul Oh, and Sabre Kais.  Implementation of quantum machine learning for electronic structure calculations of periodic systems on quantum computing devices.  preprint arXiv:2103.02037, 2021.
  * [8] Yudong Cao, Shuxian Jiang, Debbie Perouli, and Sabre Kais.  Solving set cover with pairs problem using quantum annealing.  Scientific reports, 6(1):1–15, 2016.
  * [9] S. Kais, K.B. Whaley, A.R. Dinner, and S.A. Rice.  Quantum Information and Computation for Chemistry.  Advances in Chemical Physics. Wiley, 2014.
  * [10] Katherine L Brown, Anmer Daskin, Sabre Kais, and Jonathan P Dowling.  Reducing the number of ancilla qubits and the gate count required for creating large controlled operations.  Quantum Information Processing, 14(3):891–899, 2015.
  * [11] Peter W. Shor.  Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer.  SIAM Journal on Computing, 26(5):1484–1509, Oct 1997.
  * [12] Lieven M. K. Vandersypen, Matthias Steffen, Gregory Breyta, Costantino S. Yannoni, Mark H. Sherwood, and Isaac L. Chuang.  Experimental realization of shor’s quantum factoring algorithm using nuclear magnetic resonance.  Nature, 414(6866):883–887, Dec 2001.
  * [13] Erik Lucero, R. Barends, Y. Chen, J. Kelly, M. Mariantoni, A. Megrant, P. O’Malley, D. Sank, A. Vainsencher, J. Wenner, and et al.  Computing prime factors with a josephson phase qubit quantum processor.  Nature Physics, 8(10):719–723, Aug 2012.
  * [14] Michael R. Geller and Zhongyuan Zhou.  Factoring 51 and 85 with 8 qubits.  preprint arxiv:1304.0128, 2013.
  * [15] Shuxian Jiang, Keith A. Britt, Alexander J. McCaskey, Travis S. Humble, and Sabre Kais.  Quantum annealing for prime factorization.  Scientific Reports, 8:17667, 2018.
  * [16] Nikesh S. Dattani and Nathaniel Bryans.  Quantum factorization of 56153 with only 4 qubits, 2014.
  * [17] Wilhelm Magnus.  On the exponential solution of differential equations for a linear operator.  Communications on Pure and Applied Mathematics, 7(4):649–673, 1954\.
  * [18] Yuta Shingu, Yuya Seki, Shohei Watabe, Suguru Endo, Yuichiro Matsuzaki, Shiro Kawabata, Tetsuro Nikuni, and Hideaki Hakoshima.  Boltzmann machine learning with a variational quantum algorithm, 2020\.
  * [19] Mario Motta, Chong Sun, Adrian T. K. Tan, Matthew J. O’Rourke, Erika Ye, Austin J. Minnich, Fernando G. S. L. Brandão, and Garnet Kin-Lic Chan.  Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution.  Nature Physics, 16(2):205–210, Nov 2019.
  * [20] Sam McArdle, Tyson Jones, Suguru Endo, Ying Li, Simon C. Benjamin, and Xiao Yuan.  Variational ansatz-based quantum simulation of imaginary time evolution.  npj Quantum Information, 5(1), Sep 2019.
  * [21] Kübra Yeter-Aydeniz, Bryan T. Gard, Jacek Jakowski, Swarnadeep Majumder, George S. Barron, George Siopsis, Travis S. Humble, and Raphael C. Pooser.  Benchmarking quantum chemistry computations with variational, imaginary time evolution, and Krylov space solver algorithms.  Advanced Quantum Technologies, 4:2100012.
  * [22] Chris J.C. Burges.  Factoring as optimization.  Technical Report MSR-TR-2002-83, January 2002.
  * [23] Simon J Devitt, Austin G Fowler, and Lloyd CL Hollenberg.  Investigating the practical implementation of shor’s alagorithm.  In Micro-and Nanotechnology: Materials, Processes, Packaging, and Systems II, volume 5650, pages 483–494. International Society for Optics and Photonics, 2005.
  * [24] D. PETZ and C. GHINEA.  Introduction to quantum fisher information.  Quantum Probability and Related Topics, Jan 2011.
  * [25] Héctor Abraham and AduOffei.  Qiskit: An open-source framework for quantum computing, 2019.
  * [26] Xiao Yuan, Suguru Endo, Qi Zhao, Ying Li, and Simon C. Benjamin.  Theory of variational quantum simulation.  Quantum, 3:191, Oct 2019.

##  Appendix A McLachlan’s Principle

We are borrowing the derivation from [26] and including it here for
completion.

$\delta\left\lVert(\frac{d}{d\tau}+A)\ket{\psi(\tau)}\right\rVert=0$ (20)

Here we take $A$ to be Hermitian, which is true of the Hamiltonian, and work
with the variation of the norm square

$\left\lVert(\frac{d}{d\tau}+A)\ket{\psi(\theta(\tau))}\right\rVert^{2}=((\frac{d}{d\tau}+A)\ket{\psi(\theta(\tau)))})^{*}((\frac{d}{d\tau}+A)\ket{\psi(\theta(\tau)))})$
(21)

Then, taking $\theta$ to be real, this yields

$\braket{\psi_{i}(\theta(\tau))}{\psi_{j}(\theta(\tau))}\dot{\theta_{i}}\dot{\theta_{j}}+\bra{\psi_{i}(\theta(\tau)))}A\ket{\psi(\theta(\tau)))}\dot{\theta_{i}}+\bra{\psi(\theta(\tau)))}A\ket{\psi_{i}(\theta(\tau)}\dot{\theta_{i}}+\bra{\psi(\theta(\tau)))}A^{2}\ket{\psi(\theta(\tau)))}$
(22)

where we have used the notation
$\ket{\psi_{i}(\theta)}=\frac{d\ket{\psi(\theta)}}{d\theta_{i}}$. Computing
the variation for a small time interval, we get

$(\braket{\psi_{i}(\theta(\tau))}{\psi_{j}(\theta(\tau))}\dot{\theta_{j}}+\bra{\psi_{j}(\theta(\tau)))}\frac{d\ket{\psi(\theta(\tau)}}{d\theta_{i}}\dot{\theta_{j}}+\bra{\psi(\theta(\tau)))}A\ket{\psi_{i}(\theta(\tau)}+\bra{\psi_{i}(\theta(\tau)))}A\ket{\psi(\theta(\tau)))})\delta\dot{\theta_{i}}$
(23)

Then, setting the variation to vanish for independent $\theta_{i}$, equates to
the coefficients vanishing and yields

$\Re(\braket{\psi_{i}(\theta(\tau))}{\psi_{j}(\theta(\tau))})\dot{\theta_{j}}=-\Re(\bra{\psi(\theta(\tau)))}A\ket{\psi_{i}(\theta(\tau)})$
(24)

We next substitute $A=H-\bra{\psi(\theta)}H\ket{\psi(\theta)}$ and use the
fixed normalization of the state with that
$\Re(\bra{\psi(\theta)}\ket{\psi_{i}(\theta)})=0$ to yield

$\Re(\braket{\psi_{i}(\theta(\tau))}{\psi_{j}(\theta(\tau))})\dot{\theta_{j}}=-\Re(\bra{\psi(\theta))}H\ket{\psi_{i}(\theta(\tau)})$
(25)

Finally, we get $M_{ij}\dot{\theta}=C_{i}$ with
$M_{ij}=\Re(\braket{\psi_{i}(\theta(\tau))}{\psi_{j}(\theta(\tau))})$ and
$C_{i}=\Re(\bra{\psi(\theta))}H\ket{\psi_{i}(\theta(\tau)})$

##  Appendix B Computations per iteration of the QITE

Let the number $N$ to factorized be represented by $n$ bits. Each of the
factors will be upper bounded by an $n$ bit representation. The Hamiltonian
thus takes the following form

$H=[N-(2^{n-1}p_{n-1}+..+2^{0}p_{0})(2^{n-1}q_{n-1}+..+2^{0}q_{0})]^{2}$ (26)

Expanding this Hamiltonian, we obtain a constant term, $n^{2}$ terms of the
type $p_{i}q_{j}$, $n^{2}(n-1)/2$ terms of the form $p_{i}q_{j}q_{k}$,
$n^{2}(n-1)/2$ terms of the form $p_{i}p_{j}q_{k}$ and $n^{2}(n-1)^{2}/4$
terms of the form $p_{i}p_{j}q_{k}q_{l}$. The constructed Hamiltonian thus
contains $[n(n+1)/2]^{2}+1\approx O(n^{4})$ terms in it. Note that the cost of
computing the coefficients is also of the same order, given the number of
terms is of $O(n^{4})$, if one were to ignore any simplifications that are to
be carried.

When we restrict to a $d$-layer Ry circuit, we update $nd$ circuit parameters
with every algorithm iteration. This involves computing an $M$ matrix of size
$nd$ and a C vector of size $nd$. Computing each entry in a C vector involves
working with each term of the Hamiltonian independently and summing them up to
the end, thus we have $O(n^{5}d)$ circuits to evaluate for C and
$O(n^{2}d^{2})$ circuits to evaluate for M. Taking the circuit depth to grow
slower than $n^{3}$, we get a total of $O(n^{5}d)$ circuits to compute every
iteration. To keep the overall run time of the QITE algorithm as polynomial,
it is important to achieve convergence to the correct solution in polynomial
number of steps.