# How Energy Sector Reacted to COVID-19 Pandemic? Empirical Evidence from an Emerging Market Economy

Authors: Daniel Ştefan Armeanu; Camelia Cătălina Joldeș; Ştefan Cristian Gherghina
Journal: Discrete Dynamics in Nature and Society
 (2022)

---

## Abstract

The European Union is facing the highest natural gas prices in 15 years, owing largely to an upward trend in electricity prices, which is also on an uphill curve. However, the rise in electricity and natural gas prices is a widespread phenomenon that is being felt not only in Europe but also globally, as economic activity resumes and energy consumption returns to prepandemic levels. Consequently, this paper investigates how COVID-19 influenced the Romanian energy market. To accomplish our goal, we used daily data for variables and market indices that characterize COVID-19 and the energy market from July 1 to December 21, 2021. The results of the GARCH (1, 1) model estimation show that the major performer in Romania’s energy allocation and supply market had the highest conditional variance. In addition, the ARDL model was chosen because of the variable integration mix (order 0 and 1), as well as the VAR and the Granger causality framework. The empirical results of ARDL models provide the first conclusion of the analysis, indicating that the number of short-term connections was greater than long-term connections, which is also explained by the presence of short episodes of high volatility recorded in the investigated time interval. Another conclusion drawn from this study is that COVID-19 cases registered in Europe and around the world have made a significant contribution to explaining the evolution of the energy market, owing to the large number of cases registered in these regions and the level of contagion transmitted from these markets to the energy market. Furthermore, based on the Granger causality test results, only one-way causal relationships were identified from the variables that capture the evolution of the COVID-9 pandemic to the yields of Romanian energy companies. The novelty of this article is the examination of the impact of COVID-19 on the energy market throughout the fourth wave of coronavirus using the GARCH framework, the ARDL model, which allows for the capture of both short- and long-term reactions, the variance decomposition, and the Granger causality test. Because of the ongoing changes in the pandemic’s evolution, additional research on this topic is undoubtedly on the horizon in the near future.


---

## Body

## 1. Introduction

Because of the lack of confidence, the pandemic might have an impact on international markets, causing asset prices to collapse across markets [1]. Since COVID-19 spreads quickly and has driven intensified economic insecurity [2], financial markets have plummeted and become extremely unpredictable, with significant drops in oil and metal prices [3, 4]. Szczygielski et al. [5] asserted that COVID-19 uncertainty has wreaked havoc on every national energy market. Naeem [6] supported that during times of crisis, market inefficiencies in the energy markets are more noticeable. As well, Zhang et al. [7] reinforced that the coronavirus calamity triggered considerable supply and demand shocks throughout the crude oil market. Although the oil supply and demand shocks related to COVID-19 are expected to be transitory, [8] their consequences on various sectors and nations are ongoing. In this regard, Wu and Ma [9] exhibited that COVID-19’s energy price oscillations have a harmful effect on economic growth and inflation, with a larger impact on the latter variations. Following the outbreak of coronavirus, there have been several disruptions to economic activity leading to declining production and consumption. As well, individuals were inclined to remain at home when restriction measures were adopted, causing a substantial decline in trade and manufacturing activities, which has an instant adverse effect on electricity usage [10]. Consequently, certain economic sectors were more affected, such as tourism, industry, and transport. According to Phillips [11], it is expected that the COVID-19 virus will grow endemic, but it may become less harmful over time. The return of global activity, as the restrictions imposed by the pandemic were lifted and the economies were completely reopened, led to a quick increase in demand for natural gas, both for electricity generation and for industrial purposes. In this regard, Yu et al. [12] claimed that throughout critical situations, natural resource commodities have shown to be extremely volatile. Also, the dramatic drop in crude oil prices has a substantial adverse effect on the low-carbon economy [13]. Hence, the carbon futures market was as well distorted. For instance, Dou et al. [14] emphasized that the carbon trading market has been substantially affected by the high-frequency price oscillations of carbon assets. Duan et al. [15] argued that when the economy is powerful and carbon prices are soaring, the effects of energy costs on carbon prices are less severe than when the market is fragile and carbon prices are the low point.During the summer of 2021, this condition was amplified by the strong demand for electricity production, against the background of heat waves, while reducing production from competing sources (hydro and wind) and the appearance of supply constraints, caused by extreme weather events or prolonged maintenance work (given the postponement of some of the latter in the acute phases of the pandemic). According to Li et al. [10], an unexpected variation in energy demand can have a detrimental effect on both energy services and economic decisions because it provokes issues in energy infrastructure consistency and oscillations in energy distribution systems. Specifically, in Romania, the final energy consumption was 46125.7 million kWh from January to October 2021, up 5.3 percent over the same period in 2020, while public lighting climbed by 2.6 percent and private consumption increased by 7.4 percent [16]. Bahmanyar et al. [17] reported that weekday consumption was significantly reduced in states with severe constraints (Spain, Italy, Belgium, and the United Kingdom), and energy consumption patterns were comparable to prepandemic weekend profiles for the same period in 2019. Nevertheless, the decrease in electricity usage was lower in countries with less restrictive policies. By focusing on Canadian data, Khalil and Fatmi [18] reported that the regular average in-home duration of each occupation increased by about 80% throughout the disease outbreak, causing a 29% boost in energy utilization. Similarly, Surahman et al. [19] confirmed that throughout the coronavirus period, the average annual energy utilization in Indonesia’s major cities is higher than prior to the pandemic. On the contrary, Kang et al. [20] proved in the case of South Korea that through the pandemic period, most amenities energy utilization has gradually diminished. In the same vein, Wang et al. [21] confirmed that COVID-19 has reduced China’s electricity consumption by 29%.This research explores how the COVID-19 pandemic influenced the Romanian energy market throughout the fourth wave of coronavirus. The motivation for exploring Romania is depicted by the fact that it ranked first among all European member nations and sixth internationally during the fourth COVID-19 wave. As well, according to European Centre for Disease Prevention and Control [22], in the penultimate week of December 2021, Romania has the second-lowest vaccination rate in the European Union (EU), with only 41.2 percent of its 19.3 million people fully vaccinated against COVID-19. Since 2019, Romania has switched from being an electricity exporter to a net importer, despite having historically had the third lowest rate of energy import dependency in the EU due to its natural gas and oil reserves and a large power generation sector. Also, coal, hydropower, fossil gas, nuclear energy, and wind power all contribute roughly equal amounts of capacity and power generation to Romania’s electricity mix, which is one of the most balanced in the EU [23].Our paper’s specific goals are to investigate the reaction of the energy market during the pandemic period from July to December 2021. The first goal is to examine the volatility of selected variables that characterize the energy market, namely stock market indices corresponding to the energy sector (e.g., EEX-B, EEX-P, APX, and LNGI), the Romanian energy market specific index (e.g., BET-NG), as well as companies traded on the Bucharest Stock Exchange (BSE) (e.g., SNG, EL, TGN, and TEL). The GARCH model will be considered in this regard. The second goal is to investigate long-term and short-term relationships between energy market variables and the number of new COVID-19 cases in Europe and globally. The ARDL (autoregressive distributed lag) econometric model will enable us to examine such associations [24]. Third, we will approach variance decomposition using VAR (vector autoregressive models) models to further investigate the level of explanation of COVID-19 variables on the variation of selected energy market measures. Fourth, we will investigate the type of causality that has been established between the new number of COVID-19 cases and the energy market. In this vein, after all, data series have been converted to stationary series, the Granger causality test will be applied.Crude oil prices shape macroeconomic dynamics by affecting monetary policy instruments, inflation, as well as other business activity, in addition to their impact on corporate earnings [25]. For the reason that oil price shifts have such a large impact on macroeconomic performance, exploring this effect and predicting how vulnerable economic development will be during crises is critical [12]. This research adds to the existing body of knowledge in a number of ways. First, to the best of our knowledge, the study provides the first empirical evidence for the case of Romania. Prior studies were focused on a diverse worldwide sample [4, 26], several European countries [17], China [9, 13, 21, 27–32], G7 nations [33], Indonesia [19], South Korea [20], United States [32, 34–36], Germany and United States [10], Canada [18], United Kingdom [37], United States and China [38], United States and Japan [39], China and Nigeria [40], Turkey [41], various countries with the leading energy sectors by market capitalization [5], emerging economies [2], or advanced and developing nations [42]. Hence, the evidence for emerging market economies is limited. Secondly, the most recent COVID-19 wave is covered, namely the period from July 1 to December 21, 2021. Thirdly, different from prior papers that employed merely the worldwide pandemic cases [3, 24, 43, 44], our econometric investigation covers both new cases of COVID-19 pandemic globally and in Europe. As well, contrary to prior papers that used the price of WTI crude oil [1–4, 7, 8, 24–26, 29, 33, 35, 39, 43, 45–56], Brent crude oil [4, 7, 8, 15, 24, 26, 36, 37, 43, 44, 46, 47, 49, 51, 52, 54–56], Dubai crude oil [4, 8], NYMEX’s oil [26, 55, 57], or carbon futures [14], the current study covers the daily returns of the Physical Electricity Index, Amsterdam Power Exchange Electricity Netherlands Average All Hours and, London Natural Gas Index United Kingdom Pence Per 100000 British Thermal Units. Not least, different from prior papers focused on time-frequency connectedness [32, 42, 55, 56, 58], our quantitative framework covers several techniques such as GARCH estimation, autoregressive distributed lag (ARDL) models, as well as vector autoregressive (VAR) models. Aside from empirical contributions, our findings have practical implications for policymakers throughout the unprecedented phase of insecurity triggered by the COVID-19 pandemic.This paper is organized as follows. The second section examines the previous literature. The third section presents the research sample and the selected variables, along with the applied quantitative methods. The fourth section discusses the empirical results, and the final section concludes the article.


## 2. Background Literature

Many other regions of the world have experienced a dramatic increase in energy prices during the period under review. This may be due to the increase in global energy demand (especially natural gas), as the process of economic recovery has intensified since the peak of the COVID-19 pandemic. An elementary explanation, demonstrable in numbers, is that the demand for energy in 2021 has grown too fast compared to the increase in supply. However, the explanations for this deficit are numerous: in 2020, investment and maintenance works in electricity or natural gas capacities were stopped, especially LNG in the USA; industrial demand increased significantly in early 2021 after the 2020 lockdown periods; demand in Asia has recovered from other regions and attracted a large share of global liquefied gas supplies; the spring of 2021 was colder and more gas and electricity was consumed for heating; there have been gas production problems in the US due to Hurricane Ida, and so on.It was expected that the economic activities of all households would be completely stopped throughout the world during the COVID-19 period. Most people were locked up in their homes and confined without work, which led to a loss of income. Additionally, living in homes during the pandemic has increased the electricity bill. Therefore, the situation of loss of income and increased electricity bills was a huge economic and financial burden for households [59]. Furthermore, concern for a worldwide recession has generated unavoidable systemic risks in the energy markets, exposing investors who own oil-derived securities to detrimental changes in crude prices [46]. The study by Zhang and Wang [57] demonstrated that the disease has increased long-term volatility for all future returns. Also, Shaikh [52] proved that through disease eruptions, the WTI crude oil market exhibited exceptional overreaction and dealt at an extremely volatile level. Dutta et al. [53] noticed that after the events related with COVID-19, there was a sizeable decline in worldwide crude oil prices, with the effect being highest when this novel coronavirus infection was announced a pandemic.Crude oil is sometimes considered an economy’s blood, and as a reason, oil price variations have a significant impact on many countries throughout the world [24]. For instance, Tong et al. [27] claimed that jumps in the oil market were the most strongly connected to the disease, notably through the peak and refall stages. The COVID-19 pandemic has put the energy markets under stress, but its most direct impact was on energy consumers through isolation measures that have deepened preexisting energy poverty problems, increased residential demand due to increased occupancy, and reduced the earnings of many families that have been economically affected by the crisis [60]. Narayan [61] noticed that the COVID-19 period is defined by one of the most severe drops in oil prices, with the price reaching a new low and the volatility rising by up to 900 percent.The impact on the energy sector has led to an increased demand for residential energy as a result of reduced mobility and a change in the nature of work. Blockades around the world at the beginning of the pandemic have restricted movement and placed people at home, which in turn has reduced the demand for industrial and commercial energy, as well as waste generation [62].The first strand of literature explored the impact of coronavirus diseases on energy stock returns. For example, Huang and Liu [30] found that after the pandemic period, the risk of Chinese energy companies’ stock price fall down has considerably lessened. Further, the relationship between the COVID-19 pandemic, the oil yield, and the profitability of stocks in a unified framework was achieved using a TVP-VAR model by Liu et al. [35]. There was noticed a negative association between crude oil and stock yields during the sampling period. Contrary to many people’s beliefs, the outbreak of the COVID-19 pandemic could have a significant positive impact on the crude oil market and the stock market. In this regard, Aloui et al. [63] argued that the S&P GS indexes of energy markets respond to the COVID-19 crisis in various ways throughout time due to fundamental and behavioral characteristics.Pavlyshenko [64] studied different regression approaches to model the spread of COVID-19 and its impact on the stock market. The logistic curve model was used with Bayesian regression for predictive analysis of coronavirus spread. The impact of COVID-19 was examined using regressions compared to other effects of the crisis. Empirical results showed that different crises with various causes have a distinct impact on the same stocks.In the same way, a statistical analysis of the effect of the COVID-19 pandemic on stock market risk was also conducted by Zhang et al. [65]. The virus has killed thousands of people and brought significant challenges to countries around the world. The results showed that the risks to the global financial market have increased substantially in response to the pandemic. For instance, Wang et al. [36] supported the transmission of risk among the coal and WTI crude oil markets. Individual reactions in the stock market are clearly related to the severity of the outbreak in each country. The high uncertainty of the pandemic and the associated economic losses have made markets extremely volatile and unpredictable. Political reactions are required to fight the virus and the level of stock markets. However, unconventional policy interventions, such as US quantitative easing (QE), create additional uncertainty and could cause long-term complications. Furthermore, countries do not work jointly to meet these challenges, as markets react differently to national policies and the overall development of the pandemic. Hence, this trend of fragmentation in the global community is more of a threat than a virus.The new coronavirus has generated significant volatility in financial markets, including the commodity market. Argued by the fact that oil prices have fallen the most since 1991, the second strand of literature was oriented toward the volatility investigation of oil markets throughout the pandemic period. The time-varying total volatility spillovers across markets were noticed to have strengthened with the occurrence of COVID-19 and global crude oil price turbulence throughout the pandemic [31]. Ashok et al. [26] exhibited increased co-movements in energy markets appearing months prior to co-movements in equity markets. For instance, Akyildirim et al. [42] explored 29 developed and developing states and found that oil-exporting nations predominantly spread shocks, while oil-importing states mostly receive shocks. Si et al. [28] claimed that the COVID-19 shock is one of the main factors for the Chinese energy markets to become even more volatile. Likewise, Yousaf [47] supported that the volatility of the WTI oil market increases as the volatility of COVID-19 volatility rises. Also, Wang et al. [34] found that the price of oil in the United States has fallen along with the number of new instances that have been confirmed. Albulescu [66] studied how COVID-19 figures, in terms of daily announcements of new infections, have influenced international oil prices. The ARDL estimate showed a negative and significant impact on the coronavirus crisis but was relatively small compared to the effect of financial volatility and uncertainty in economic policy on oil prices. However, the influence of the pandemic on oil prices was indirect, and the volatility of mainly affected financial market. In addition, Albulescu [67] emphasized that the downward fluctuations of crude prices are driven by increased insecurity. The outcomes were also supported by Jeris and Nath [37], and Geyikçi [41]. Narayan [61] reinforced that the oil market has become unprofitable over the COVID-19 time frame using technical moving average trading techniques. Lin and Su [50] found that subsequent to the eruption of COVID-19, there is a dramatic rise in total connectedness in energy markets, but this shift merely persists about two months prior to returning to preinfection levels. Therefore, Iglesias and Rivera-Alonso [51] argued that volatility peaks occur throughout periods of supply/demand downturns or oil instabilities, whereas cycles with financial turmoil as the main trigger are associated with higher volatility persistence.Other studies explored the impact of pandemic news on the energy sector. For instance, Albulescu [68] examined the impact of official COVID-19 announcements and related statistics on financial volatility, comparing the effect of data reported in China with those of COVID-19 records reported outside China. Empirical results exposed that only new cases reported outside China have a positive effect on the VIX index. Also, the death rate has a positive influence on the VIX index for all estimated models, but the effect was greater for the death rate outside China. In addition, the spread of the pandemic increased financial volatility. As such, the persistence of COVID-19 could generate a new episode of global financial stress. Amamou and Bargaoui [49] found that the release of a new disease outbreak wave lessens dependence on the oil market, which losses its attributes as a safe-haven market in favor of other markets such as gold or cryptocurrencies. Shaikh [54] exhibited that global crude oil is negatively associated with the news connected to the COVID-19 pandemic.Further, Akhtaruzzaman et al. [69] investigated how financial contagion occurs through financial and nonfinancial firms between China and the G7 countries. Empirical results revealed that the dynamic conditional correlation (DCC) between the profitability of Chinese financial stocks and the financial and nonfinancial G7 increased significantly during the COVID-19 period. However, the magnitude of DCC growth has been greater for financial firms, implying that they exert a more critical role in transmitting the financial contagion than nonfinancial firms. The results showed that China and Japan appear to be net emitters of COVID-19 contagion. In addition, optimal hedging ratios increased substantially in most cases during the COVID-19 phase, leading to higher hedging costs during the crisis. The findings of Jiang and Chen [55] also indicate that overall connectedness has increased significantly since the COVID-19 outbreak, as compared to the preCOVID era. Jiang and Chen [58] confirmed that following the COVID-19 outbreak, overall connectivity has increased. Equally, Mensi et al. [38] exhibited that during the low-volatility period (high volatility regime), oil was a prominent recipient (contributor) of spillovers. Chen et al. [32] remarked that prior to the COVID-19 pandemic, the role of energy commodities as net receivers can be noticed in both the short and long terms, but following it, the net transmitter position can be seen in the long term.Further, a review of earlier research exploring the influence of the COVID-19 pandemic on energy markets is provided in Table1.Table 1
Summary of prior literature on the COVID-19–energy market nexus.
Author (s)PeriodVariablesMethodologyOutcomesMa et al. [29]January 1, 2019–April 1, 2021Oil prices and the GDPWavelet power spectrum, wavelet coherence, frequency domain causality testThe price of natural resource commodities has been observed to be more volatile throughout the COVID-19 timeframeYu et al. [12]2007–2009, 2010–2018, 2019–2020GDP and oil pricesWavelet analysisIn the long-term, a jump in crude prices has effects on economic growthAli et al. [25]March 2020–May 2020Closing spot prices of WTI crude oil futures and stock indices of the United States, Canada, China, Russia, and VenezuelaWavelet-based granger causalityThroughout times of stability, oil is vital for hedging, and during times of crisis, it serves as a safe-haven assetMensi et al. [70]April 23, 2018–April 24, 2020S&P500 index, brent oil, and gold futuresBivariate FIAPARCH modelFor all sub-periods, oil offers greater hedging efficiency than goldAtri et al. [3]January 23, 2020–June 23, 2020WTI oil price, the worldwide confirmed new cases, and deathsARDL analysisDuring the COVID-19 contagion, economic and financial instability has a detrimental effect on oil and gold valuesGharib et al. [71]November 1, 2019–December 31, 2020Daily west Texas light crude oil and north sea brent crude, diesel, and gasoline pricesSupremum augmented Dickey-Fuller, generalized supremum augmented Dickey-Fuller, the explosive test strategyThroughout the COVID-19 outbreak, west Texas Light crude oil and north sea brent crude oil had a negative financial bubbleAdedeji et al. [40]March 20, 2020–May 28, 2020West Texas intermediate, brent, bonny, DaqingVector autoregressive (VAR) methodThe influence of the COVID-19 pandemic on bonny and Daqing oil prices accounted for the smallest shares of fluctuation, while the effect on BRENT and WTI is even smallerBourghelle et al. [48]January 2, 2014–April 1, 2020West Texas index, economic policy uncertainty index, equity market-related EPU indexVector autoregressive (VAR) frameworkThe disease outbreak oil shocks had a significant effect on oil price fluctuationsKhan et al. [72]January 2020–May 2021West intermediate Texas, brent oil, natural gas, heating oilQuantile-on-quantile methodCOVID-19 has a generally negative impact on energy prices across all quantilesManeejuk et al. [44]December 29, 2019–December 30, 2020Natural gas, gasoline, heating oil, coal, and brent crude oilGeneralized autoregressive conditional heteroskedasticity (GARCH), Markov switching dynamic copulaEnergy markets react the same way to both positive and negative occurrences of COVID-19Nyga-Łukaszewska and Aruga [39]January 1, 2020–June 2, 2020,WTI, Platts Dubai crude oil prices, Henry hub, Platts Japan Korea marker pricesAuto-regressive distributive lag (ARDL)The COVID-19 pandemic in the United States had a statistically detrimental effect on crude oil prices while having a positive impact on gas pricesLe et al. [24]January 17, 2020–September 14, 2020WTI oil price, brent, trade-weighted US dollar index, MSCI world index, FTSE all-world index, S&P Global 100 indexARDL bounds testing procedureThe fall in WTI prices is attributed to increases in COVID-19 instances, US economic policy uncertainty index, and the Chicago board options exchange (CBOE) volatility index (VIX)Ahundjanov et al. [43]January 22, 2020–July 2, 2020Brent, west Texas intermediate, New York harbor, Dow Jones US oil and gasStructural vector autoregressive (SVAR) modelA unit rise in COVID-19 global search interest leads to a cumulative reduction of 0.083 percent and 0.104 percent in the Dow Jones US oil, and gas total index and New York harbor conventional gasoline, respectivelyLi et al. [33]December 1, 2019–March 25, 2022.WTI crude oil futures pricesMultivariate waveletThe COVID-19 pandemic is less of a concern to the people in the US and Canada than the fall in the WTI and worldwide stock marketsChatziantoniou et al. [56]January 17, 1997–December 11, 2020WTI, brent, heating oil, kerosene, propane, and gasolineConditional autoregressive value-at-risk (CAViaR), TVP-VARWith substantial crisis events, connectivity increasesSource: author’s own work.


## 3. Research Methodology

### 3.1. Sample and Variables

The European Union, similar to various other regions of the world, is currently facing a sharp rise in energy prices. The world economy is significantly impacted by the pandemic caused by COVID-19 disease, and the effects will be long-term and will depend on the intensity of the pandemic. In our quantitative study, we use daily data from July 1, 2021, to December 21, 2021, namely the period that incorporates the fourth wave of the COVID-19 pandemic.To achieve the proposed goals, we picked the number of new COVID-19 cases registered in Europe and globally as COVID-19 variables. We assume that the evolution of the number of new COVID-19 cases will better capture the pandemic’s impact on the energy market, where we used various indices that seize the price of electricity and gas. In addition, we intended to include a Romanian stock market index that tracks the evolution of companies in the energy sector, as well as a survey of the composite index’s largest companies. The selected measures are presented in Table2.Table 2
Variables’ description.
VariablesDescriptionCNWNew cases of COVID-19 pandemic worldwideCNENew cases of COVID-19 pandemic in EuropeBET-NGThe daily return of BET-NG–is a sectoral index that reflects the evolution of companies listed on the regulated market of the Bucharest stock exchange that have the main field of activity energy and related utilitiesSNGThe daily return of S.N.G.N. ROMGAZ S.A.–company that has a vast experience in the field of exploration and production of natural gas, being one of the largest producers and main suppliers of natural gas in RomaniaELThe daily return of electrica–is the only company in Romania listed in the field of electricity distribution and supplyTGNThe daily return of SNTGN Transgaz SA–is the technical operator of the national transport system (NTS) natural gas that ensures conditions of safety, efficiency, competitiveness, and in compliance with European standards of performance and environment, the transport of over 90% of natural gas consumed in RomaniaTELThe daily return of CNTEE trans electrica SA–is already recognized on a national scale and globally as a strong company, with a strategic role in the Romanian electricity market and an essential participant in the regional electricity marketEEX-BThe daily return of the physical electricity index (Phelix)–refers to the base load (Phelix base) price index published daily on the power spot market for the German/Austrian market areaEEX-PThe daily return of the physical electricity index (Phelix)–refers to the peak load (Phelix peak) price index published daily on the power spot market for the German/Austrian market areaAPXThe daily return of the Amsterdam power exchange (APX) electricity Netherlands average all hoursLNGIThe daily return of the London natural gas index United Kingdom pence per 100000 British thermal unitsSource: author’s own work.Therefore, a wide array of variables has been selected from the Refinitiv Eikon database alike Ashok et al. [26], which allows the achievement of our research objective, namely, indices on the stock market that describe the energy sector (e.g., EEX-B, EEX-P, APX, and LNGI), along with the specific Romanian energy market index BET-NG. In addition, we have included the daily returns of several leading Romanian companies listed on the Bucharest Stock Exchange (BSE) that operate in the energy field (e.g., SNG, EL, TGN, and TEL). COVID-19 pandemic is measured through daily recent reported cases worldwide (CNW) alike Villarreal-Samaniego [2], Atri et al. [3], Le et al. [24], Jeris and Nath [37], Ahundjanov et al. [43], Maneejuk et al. [44] and in Europe (CNE).


### 3.2. Quantitative Analysis Strategy

To explore how the COVID-19 pandemic affects the Romanian energy market, more specifically the price of energy and the price of gas, we will employ a variety of econometric tools such as: (1) stationarity analysis, (2) generalized autoregressive conditional heteroskedasticity (GARCH) estimation, (3) autoregressive distributed lag (ARDL) models, (4) vector autoregressive (VAR) models, (5) the Granger causality test.First, the Augmented Dickey–Fuller Unit Root (ADF) test will be used to verify the nonstationarity of our variables similar to Villarreal-Samaniego [2], Atri et al. [3], Bildirici et al. [8], Wu and Ma [9], Li et al. [10], Wang et al. [21], Wang et al. [34], Geyikçi [41], Ahundjanov et al. [43], Amamou and Bargaoui [49], Lin and Su [50], Zhang and Wang [57], Albulescu [66], Albulescu [67]. Specifically, nonstationary variables lead to inadequate results, which means insignificant results. The confirmation of the stationarity of the selected data is performed through the ADF stationarity test, being the regular test employed to confirm the stationarity of a data series. The null hypothesis of the ADF test assumes that the variable has a unit root, and thus the measure is not stationary. The ADF test involves estimating the equation as follows:(1)Δmt=α+βt+qωt+∑j=1kγjΔmt−j+εt,t=1,…,T,where t represents the time trend, T is sample length and k is the length of the lag in the dependent variable. Nevertheless, the ADF test is based on a linear assumption, which can lead to inaccurate results [8]. Thus, to check the robustness of our results, further unit root tests alike Kwiatkowski-Phillips-Schmidt-Shin and Zivot–Andrews with one break will be applied in line with Zhang, Farnoosh [7].Secondly, volatility clustering is a method of identifying market volatility. Due to the pandemic news outburst, the market typically encounters an unstable phase prior to returning to regularity [52]. In order to explore the volatility of the selected series, the GARCH (p,q) model will be considered alike [4], Szczygielski et al. [5], Wu and Ma [9], Maneejuk et al. [44], Iglesias and Rivera-Alonso [51], Zhang and Wang [57], which shows the following general form:(2)σt2=α0+∑i=1qαiut−i2+∑j=1pβjσt−j2,where p is the order of the GARCH terms and q is the order of the ARCH terms.Thirdly, in case we find out both stationary and nonstationary variables, consistent with Villarreal-Samaniego [2], Atri et al. [3], Li et al. [10], Le et al. [24], Jeris and Nath [37], Nyga-Łukaszewska and Aruga [39], Geyikçi [41], Albulescu [66], Albulescu [67], the ARDL approach will allow the study of both short-term and long-term relationships between COVID-19 variables and the energy market. The ARDL approach exhibits several advantages over different cointegration models in quantitative literature [37]. First, the autoregressive distributed lag model (ARDL) and the limit testing methodology will be used due to its permission to apply a mixture of variables I (0) and I (1) [2, 37, 41]. Choosing the appropriate ARDL model will allow us exploring the relationships that are established between variables, so it is imperative to select the proper number of offsets. Therefore, the Akaike information criteria (AIC) will be examined to decide on the optimal gaps for the variables included in the ARDL model alike [2]. Second, this method abridges the study of the link among the response and input variables through OLS regressions. Third, when contrasted to other methodologies, the ARDL technique is more effective for small samples, which is particularly critical for this research [2, 37, 41]. Not least, the ARDL specification permits simultaneous estimation of both long-run and short-run parameters [37, 41].Specifically, anARDL (p, q1, …, qk) is a least squares regression containing lags of the dependent (p) and explanatory variables (q1, …, qk). The general specification of an ARDL (p,q) model is depicted below:(3)Ht=μ+β0Kt+β1Kt−1+⋯+βqKt−q+δ1Ht−1+⋯+δpHt−p+ut.Fourthly, the vector autoregression (VAR) framework is considered due to its common practice for interdependent time series prediction systems and for analyzing the dynamic impact of random perturbations on the system of variables. The equation for the VAR model is depicted below, consistent with Wu and Ma [9], Bourghelle et al. [48]:(4)xt=a1xt−1+⋯+apxt−p+bnt+ϵt,where xt is a k vector of endogenous variables, nt is a d vector of exogenous variables, a1, …, ap and b are matrices of coefficients to be estimated, and ϵt is a vector of innovations.Not least, causality between variables will be examined through the Granger causality test, like Wu and Ma [9], Bourghelle et al. [48]. The null hypothesis of the test consists of the following statements: h does not cause Granger k and that k does not cause Granger h. There are estimated the following bivariate regressions:(5)kt=α0+α1kt−1+⋯+αpkt−p+β1ht−1+⋯+βph−p+ϵt,ht=α0+α1ht−1+⋯+αpht−p+β1kt−1+⋯+βpk−p+ut.


## 3.1. Sample and Variables

The European Union, similar to various other regions of the world, is currently facing a sharp rise in energy prices. The world economy is significantly impacted by the pandemic caused by COVID-19 disease, and the effects will be long-term and will depend on the intensity of the pandemic. In our quantitative study, we use daily data from July 1, 2021, to December 21, 2021, namely the period that incorporates the fourth wave of the COVID-19 pandemic.To achieve the proposed goals, we picked the number of new COVID-19 cases registered in Europe and globally as COVID-19 variables. We assume that the evolution of the number of new COVID-19 cases will better capture the pandemic’s impact on the energy market, where we used various indices that seize the price of electricity and gas. In addition, we intended to include a Romanian stock market index that tracks the evolution of companies in the energy sector, as well as a survey of the composite index’s largest companies. The selected measures are presented in Table2.Table 2
Variables’ description.
VariablesDescriptionCNWNew cases of COVID-19 pandemic worldwideCNENew cases of COVID-19 pandemic in EuropeBET-NGThe daily return of BET-NG–is a sectoral index that reflects the evolution of companies listed on the regulated market of the Bucharest stock exchange that have the main field of activity energy and related utilitiesSNGThe daily return of S.N.G.N. ROMGAZ S.A.–company that has a vast experience in the field of exploration and production of natural gas, being one of the largest producers and main suppliers of natural gas in RomaniaELThe daily return of electrica–is the only company in Romania listed in the field of electricity distribution and supplyTGNThe daily return of SNTGN Transgaz SA–is the technical operator of the national transport system (NTS) natural gas that ensures conditions of safety, efficiency, competitiveness, and in compliance with European standards of performance and environment, the transport of over 90% of natural gas consumed in RomaniaTELThe daily return of CNTEE trans electrica SA–is already recognized on a national scale and globally as a strong company, with a strategic role in the Romanian electricity market and an essential participant in the regional electricity marketEEX-BThe daily return of the physical electricity index (Phelix)–refers to the base load (Phelix base) price index published daily on the power spot market for the German/Austrian market areaEEX-PThe daily return of the physical electricity index (Phelix)–refers to the peak load (Phelix peak) price index published daily on the power spot market for the German/Austrian market areaAPXThe daily return of the Amsterdam power exchange (APX) electricity Netherlands average all hoursLNGIThe daily return of the London natural gas index United Kingdom pence per 100000 British thermal unitsSource: author’s own work.Therefore, a wide array of variables has been selected from the Refinitiv Eikon database alike Ashok et al. [26], which allows the achievement of our research objective, namely, indices on the stock market that describe the energy sector (e.g., EEX-B, EEX-P, APX, and LNGI), along with the specific Romanian energy market index BET-NG. In addition, we have included the daily returns of several leading Romanian companies listed on the Bucharest Stock Exchange (BSE) that operate in the energy field (e.g., SNG, EL, TGN, and TEL). COVID-19 pandemic is measured through daily recent reported cases worldwide (CNW) alike Villarreal-Samaniego [2], Atri et al. [3], Le et al. [24], Jeris and Nath [37], Ahundjanov et al. [43], Maneejuk et al. [44] and in Europe (CNE).


## 3.2. Quantitative Analysis Strategy

To explore how the COVID-19 pandemic affects the Romanian energy market, more specifically the price of energy and the price of gas, we will employ a variety of econometric tools such as: (1) stationarity analysis, (2) generalized autoregressive conditional heteroskedasticity (GARCH) estimation, (3) autoregressive distributed lag (ARDL) models, (4) vector autoregressive (VAR) models, (5) the Granger causality test.First, the Augmented Dickey–Fuller Unit Root (ADF) test will be used to verify the nonstationarity of our variables similar to Villarreal-Samaniego [2], Atri et al. [3], Bildirici et al. [8], Wu and Ma [9], Li et al. [10], Wang et al. [21], Wang et al. [34], Geyikçi [41], Ahundjanov et al. [43], Amamou and Bargaoui [49], Lin and Su [50], Zhang and Wang [57], Albulescu [66], Albulescu [67]. Specifically, nonstationary variables lead to inadequate results, which means insignificant results. The confirmation of the stationarity of the selected data is performed through the ADF stationarity test, being the regular test employed to confirm the stationarity of a data series. The null hypothesis of the ADF test assumes that the variable has a unit root, and thus the measure is not stationary. The ADF test involves estimating the equation as follows:(1)Δmt=α+βt+qωt+∑j=1kγjΔmt−j+εt,t=1,…,T,where t represents the time trend, T is sample length and k is the length of the lag in the dependent variable. Nevertheless, the ADF test is based on a linear assumption, which can lead to inaccurate results [8]. Thus, to check the robustness of our results, further unit root tests alike Kwiatkowski-Phillips-Schmidt-Shin and Zivot–Andrews with one break will be applied in line with Zhang, Farnoosh [7].Secondly, volatility clustering is a method of identifying market volatility. Due to the pandemic news outburst, the market typically encounters an unstable phase prior to returning to regularity [52]. In order to explore the volatility of the selected series, the GARCH (p,q) model will be considered alike [4], Szczygielski et al. [5], Wu and Ma [9], Maneejuk et al. [44], Iglesias and Rivera-Alonso [51], Zhang and Wang [57], which shows the following general form:(2)σt2=α0+∑i=1qαiut−i2+∑j=1pβjσt−j2,where p is the order of the GARCH terms and q is the order of the ARCH terms.Thirdly, in case we find out both stationary and nonstationary variables, consistent with Villarreal-Samaniego [2], Atri et al. [3], Li et al. [10], Le et al. [24], Jeris and Nath [37], Nyga-Łukaszewska and Aruga [39], Geyikçi [41], Albulescu [66], Albulescu [67], the ARDL approach will allow the study of both short-term and long-term relationships between COVID-19 variables and the energy market. The ARDL approach exhibits several advantages over different cointegration models in quantitative literature [37]. First, the autoregressive distributed lag model (ARDL) and the limit testing methodology will be used due to its permission to apply a mixture of variables I (0) and I (1) [2, 37, 41]. Choosing the appropriate ARDL model will allow us exploring the relationships that are established between variables, so it is imperative to select the proper number of offsets. Therefore, the Akaike information criteria (AIC) will be examined to decide on the optimal gaps for the variables included in the ARDL model alike [2]. Second, this method abridges the study of the link among the response and input variables through OLS regressions. Third, when contrasted to other methodologies, the ARDL technique is more effective for small samples, which is particularly critical for this research [2, 37, 41]. Not least, the ARDL specification permits simultaneous estimation of both long-run and short-run parameters [37, 41].Specifically, anARDL (p, q1, …, qk) is a least squares regression containing lags of the dependent (p) and explanatory variables (q1, …, qk). The general specification of an ARDL (p,q) model is depicted below:(3)Ht=μ+β0Kt+β1Kt−1+⋯+βqKt−q+δ1Ht−1+⋯+δpHt−p+ut.Fourthly, the vector autoregression (VAR) framework is considered due to its common practice for interdependent time series prediction systems and for analyzing the dynamic impact of random perturbations on the system of variables. The equation for the VAR model is depicted below, consistent with Wu and Ma [9], Bourghelle et al. [48]:(4)xt=a1xt−1+⋯+apxt−p+bnt+ϵt,where xt is a k vector of endogenous variables, nt is a d vector of exogenous variables, a1, …, ap and b are matrices of coefficients to be estimated, and ϵt is a vector of innovations.Not least, causality between variables will be examined through the Granger causality test, like Wu and Ma [9], Bourghelle et al. [48]. The null hypothesis of the test consists of the following statements: h does not cause Granger k and that k does not cause Granger h. There are estimated the following bivariate regressions:(5)kt=α0+α1kt−1+⋯+αpkt−p+β1ht−1+⋯+βph−p+ϵt,ht=α0+α1ht−1+⋯+αpht−p+β1kt−1+⋯+βpk−p+ut.


## 4. Empirical Findings

### 4.1. Summary Statistics

The descriptive statistics of the variables are provided in Table3. The skewness and kurtosis indicators indicate the deviation in relation to a symmetric distribution around the average, thus suggesting the degree of flattening or sharpening. A kurtosis greater than the value of three implies that the returns of the indices show heavy tails than the normal distribution. Specifically, the probability of extreme returns is higher than the probability that they are below normal distribution. This feature is called leptokurtic or basically heavy tails. A positive skewness signifies an asymmetric distribution on the right and a negative skewness on the left. For a series with a normal distribution, kurtosis takes the value of three. However, for a value less than three the distribution is flatter than the normal one (e.g., platykurtic), whereas for kurtosis greater than three, the distribution is leptokurtic. According to Table 3, most of the variables exhibit a value of the kurtosis greater than three, thus presenting a leptokurtic distribution, consistent with Bildirici et al. [8], Wu and Ma [9], Wang et al. [36], Geyikçi [41], Bourghelle et al. [48], Dutta et al. [53]. Consequently, extreme negative returns are much more likely to occur than normal distribution forecasts.Table 3
Descriptive statistics of the variables.
VariablesMeanStd. dev.SkewnessKurtosisJarque-BeraProbabilityCNE2310541089950.892.4418.010.00CNW5788641130850.352.384.500.11BET-NG0.00010.01−0.744.4021.550.00SNG0.00120.01−0.826.5679.630.00EL−0.00200.01−0.755.2337.410.00TEL−0.00120.020.488.58165.980.00TGN−0.00220.01−0.374.5214.750.00EEX-B0.01240.270.155.4932.580.00EEX-P0.01340.340.105.8843.040.00APX0.01240.190.596.5371.430.00LNGI0.00990.062.0643.218442.140.00Source: author’s own work. Notes: for the definition of variables, please see Table2.The normality of the variable distribution is provided by the Jarque-Bera test. Table3 also shows the results of the Jarque-Bera test, which indicate that the distribution of the variables is not distributed normally, in line with Villarreal-Samaniego [2], Bildirici et al. [8], Ali et al. [25], Si et al. [28], Wang et al. [36], Mensi et al. [38], Akhtaruzzaman et al. [69], Gharib et al. [71], Khan et al. [72]. Because the Jarque-Bera statistic is significant (except for the number of new cases of COVID-19 registered globally), we reject the null hypothesis of normality. The test values are quite far from the corresponding normal distribution, which supports that the series is not normally distributed. Figure 1 exhibits the density, distribution, and quantile-quantile (QQ) plots. In line with the outcomes reported in Table 3, the distribution of the selected series is dissimilar from the normal one.Figure 1
Density, distribution, and QQ plots for daily returns. Source: author’s own work. Notes: variables’ descriptions are provided in Table2.The correlations between included variables are plotted in Figure2. Alike [37], the correlations were reduced. Therefore, our empirical outcomes will not be affected by the multicollinearity issue.Figure 2
Correlations among selected variables. Source: author’s own work. Notes: variables’ descriptions are provided in Table2.Figure3 shows the evolution of the number of new cases at the European and global levels due to COVID-19. Cássaro and Pires [73] argued that the number of cases is growing rapidly, which has been achieved by stability later, this mode is called a step function.Figure 3
Evolution of COVID-19 pandemic variables. Source: author’s own work. Notes: for the definition of variables, please see Table2.Further, the daily evolution of the BET-NG index and the companies listed on BSE are presented in Figure4. During the third and fourth quarters of 2021, BET-NG, the index of the 10 Romanian companies in the energy and utility sectors, registered a relatively steady evolution. EL, SNG, and TGN, which are among the most traded companies on BSE, had a similar evolution to the BET-NG index, but TEL was marked by slightly more significant episodes of volatility than the rest of the energy companies.Figure 4
Evolution of BET-NG index and of the companies listed on BSE. Source: author’s own work. Notes: for the definition of variables, please see Table2.Likewise, the evolution related to the return of the Physical Electricity Index (Phelix), Amsterdam Power Exchange (APX), and London Natural Gas Index is exhibited in Figure5. Energy prices have witnessed high episodes of volatility compared to the price of natural gas. The return of global activity, as the restrictions imposed by the pandemic and the complete reopening of economies have lifted, has led to a rapid increase in demand for natural gas, both for electricity production and for manufacturing reasons.Figure 5
Evolution of energy market. Source: author’s own work. Notes: for the definition of variables, please see Table2.


### 4.2. The Outcomes of Time Series Investigation

#### 4.2.1. Stationarity Analysis

Tables4 and 5 show the outcomes of stationarity analysis for series at the level and in the first difference through the augmented Dickey–Fuller and Kwiatkowski–Phillips–Schmidt–Shin unit root test. A couple of variables are stationary at the first difference (e.g., new cases of COVID-19 worldwide and in Europe) because the probability is above the 1% and 5% relevance level, while some are already stationary (e.g., BET-NG, SNG, EL, TGN, TEL, EEX-B, EEX-P, APX, and LNGI). Alike Bildirici et al. [8], Ali et al. [25], Si et al. [28], Mensi et al. [38], and Yousaf [47], we notice a common integration order of I (0), exception makes the indicators related to the evolution of COVID-19, so in this case, we can reject the null hypothesis and conclude that the series is not stationary.Table 4
The outcomes of the augmented Dickey–Fuller unit root test.
VariablesLevel1st differencet-statisticProb.t-statisticProb.BET-NG−10.968150.00000−8.025040.00000SNG−9.797470.00000−7.022060.00000EL−10.892070.00000−8.212170.00000TGN−11.573740.00000−8.255440.00000TEL−10.629240.00000−7.738480.00000EEX-B−9.284200.00000−10.038590.00000EEX-P−9.262510.00000−7.880230.00000APX−3.852240.00330−9.075060.00000LNGI−10.804120.00000−6.786370.00000D (CNW)−1.054550.73150−6.464690.00000CNE1.492120.99920−3.174390.02410Source: author’s own work. Notes: null hypothesis: each series has a unit root. Intercept included in test equation. Lag length: automatic selection based on Akaike info criterion. Test critical values: 1% level: −3.484198; 5% level: −2.885051; 10% level: −2.579386. For the definition of variables, please see Table2.Table 5
The outcomes of the Kwiatkowski–Phillips–Schmidt–Shin unit root test.
VariablesLevel1st differenceLM-stat.LM-stat.BET-NG0.123170.10802SNG0.051320.27447EL0.073580.22224TGN0.335890.09240TEL0.095700.03197EEX-B0.264770.10671EEX-P0.254600.10521APX0.300140.14376LNGI0.068380.33024D (CNW)0.119720.09233CNE1.172860.13543Source: author’s own work. Notes: null hypothesis: each series is stationary. Intercept included in test equation. Asymptotic critical values: 1% level: 0.739000; 5% level: 0.463000; 10% level: 0.347000. For the definition of variables, please see Table2.Further, Table6 reveals the outcomes of stationarity analysis of the variables with one structural break by means of the Zivot-Andrews unit root test, alike Zhang et al. [7], Le et al. [24]. Accordingly, the occurrence of a structural break in our sample is proved, while the mix integration is strengthened.Table 6
The outcomes of the Zivot–Andrews unit root test.
VariablesLevel1st differencet-statisticProb.Chosen breakpointt-statisticProb.Chosen breakpointBET-NG−11.602500.0145810/29/2021−10.169530.369439/28/2021SNG−10.093280.1099810/28/2021−8.848800.3967611/11/2021EL−11.146200.0322710/28/2021−7.832430.3128411/22/2021TGN−12.346240.0098311/05/2021−9.038700.1122611/24/2021TEL−10.828270.0950610/26/2021−8.099400.0481410/22/2021EEX-B−9.512800.0407710/25/2021−9.625170.1805010/22/2021EEX-P−9.478390.2147810/11/2021−11.075400.1708110/22/2021APX−8.482670.0356010/14/2021−10.477770.0891710/08/2021LNGI−11.552590.0070910/07/2021−8.209490.0262510/04/2021D (CNW)−9.573620.000028/23/2021−15.985630.342039/16/2021CNE−2.455740.0013010/19/2021−6.629530.0239310/06/2021Source: author’s own work. Notes: null hypothesis: each series has a unit root with a structural break in the intercept. Break included in the intercept. 1% critical value: −5.34.5% critical value: −4.93.10% critical value: −4.58. For the definition of variables, please see Table2.


#### 4.2.2. Volatility Examination

Because all Jarque–Bera figures are larger than three as Bildirici et al. [8], Zhang and Wang [57] found, and the distributions of the selected series diverge from normally distributed data; GARCH models are required and appropriate for highlighting the progress of volatility. Further, the autocorrelation function (ACF) and partial autocorrelation function (PACF) are used to test for heteroskedasticity. In this regard, Figure 6 exhibits the plots of ACF and PACF. Thus, the existence of serial correlation is inherent in the correlogram of squared returns, except LNGI for which the GARCH model will not be estimated.Figure 6
Plots of autocorrelation function (ACF) and partial autocorrelation function (PACF) for squared returns. Source: author’s own work. Notes: the lag length was set to 20. Variables’ descriptions are provided in Table2.Table7 reports the findings after estimating the GARCH (1, 1) model. In the case of TEL, the sum of ARCH and GARCH parameters is roughly equal to one (e.g., 0.989732), in line with Napon and Asama [4], Bildirici et al. [8]. Hence, at time t, a shock wave will last longer. Specifically, the conditional variance is persistent consistent with Maneejuk et al. [44], Iglesias, and Rivera-Alonso [51]. The findings, however, show the existence of a mean-reverting process because the sum of the ARCH and GARCH effects is less than one. However, in the case of EEX-B, EEX-P, and APX, the sum of parameters is larger than one, thus suggesting that the conditional variance process is explosive.Table 7
The outcomes of the GARCH (1, 1) model.
Dependent variable: BET-NGDependent variable: SNGVariableCoefficientStd. errorz-statisticProb.VariableCoefficientStd. errorz-statisticProb.C0.0005120.0007770.6589560.5099C0.0022260.0009642.307890.021Variance equationVariance equationC4.21E − 054.80E − 050.8768020.3806C9.67E − 050.0001810.5333310.5938RESID (−1)2−0.074930.068867−1.0880460.2766RESID (−1)2−0.048840.072398−0.6746020.4999GARCH (−1)0.6563340.4563521.438220.1504GARCH (−1)0.5420730.9169320.5911810.5544T-DIST. DOF5.0615213.1427831.6105220.1073T-DIST. DOF3.4526231.2507472.7604480.0058ARCH + GARCH0.581404ARCH + GARCH0.493233R-squared−0.001889Mean dependent var8.42E − 05R-squared−0.006586Mean dependent var0.001153AdjustedR-squared−0.001889S.D. dependent var0.009884AdjustedR-squared−0.006586S.D. dependent var0.013269S.E. of regression0.009894Akaike info criterion−6.409628S.E. of regression0.013313Akaike info criterion−5.904976Sum squared resid0.01204Schwarz criterion−6.295907Sum squared resid0.0218Schwarz criterion−5.791255Log likelihood402.3969Hannan–Quinn criterion.−6.363432Log likelihood371.1085Hannan–Quinn criterion.−5.85878Durbin-Watson stat1.988924Durbin-Watson stat1.7568Dependent variable: ELDependent variable: TGNVariableCoefficientStd. errorz-statisticProb.VariableCoefficientStd. errorz-statisticProb.C−8.80E − 060.000686−0.0128180.9898C−0.0015810.000705−2.2426210.0249Variance equationVariance equationC0.0003490.0082760.0421860.9664C6.14E − 065.86E − 061.0488260.2943RESID (−1)2−1.78439542.03255−0.0424530.9661RESID (−1)20.0757880.0746651.0150350.3101GARCH (−1)0.9877170.01889952.262890GARCH (−1)0.8743390.1022928.5474640T-DIST. DOF2.0243890.584223.4651160.0005T-DIST. DOF4.3495662.4443851.7794120.0752ARCH + GARCH−0.796678ARCH + GARCH0.950127R-squared−0.025806Mean dependent var−0.001991R-squared−0.004356Mean dependent var−0.002206AdjustedR-squared−0.025806S.D. dependent var0.012388AdjustedR-squared−0.004356S.D. dependent var0.009498S.E. of regression0.012547Akaike info criterion−6.171115S.E. of regression0.009519Akaike info criterion−6.559459Sum squared resid0.019363Schwarz criterion−6.057394Sum squared resid0.011145Schwarz criterion−6.445738Log likelihood387.6091Hannan–Quinn criterion.−6.124919Log likelihood411.6864Hannan–Quinn criterion.−6.513262Durbin–Watson stat1.93045Durbin–Watson stat2.091561Dependent variable: TELDependent variable: EEX-BVariableCoefficientStd. errorz-statisticProb.VariableCoefficientStd. errorz-statisticProb.C−0.0008560.000926−0.9240740.3554C0.0061580.0107720.5716660.5675Variance equationVariance equationC5.27E − 061.80E − 062.9369290.0033C0.0060220.0147580.4080330.6832RESID (−1)2−0.050780.016305−3.1144110.0018RESID (−1)21.519244.012510.3786260.705GARCH (−1)1.0405120.0139374.696550GARCH (−1)0.671250.1101996.0912320T-DIST. DOF5.7675783.2191841.7916270.0732T-DIST. DOF2.2005450.5998483.6685030.0002ARCH + GARCH0.989732ARCH + GARCH2.19049R-squared−0.000593Mean dependent var−0.001222R-squared−0.000556Mean dependent var0.012439AdjustedR-squared−0.000593S.D. dependent var0.015068AdjustedR-squared−0.000556S.D. dependent var0.267469S.E. of regression0.015073Akaike info criterion−5.758566S.E. of regression0.267543Akaike info criterion−0.16669Sum squared resid0.027944Schwarz criterion−5.644845Sum squared resid8.804247Schwarz criterion−0.052969Log likelihood362.0311Hannan–Quinn criterion.−5.71237Log likelihood15.3348Hannan–Quinn criterion.−0.120494Durbin-Watson stat1.930154Durbin-Watson stat2.342452Dependent variable: EEX-PDependent variable: APXVariableCoefficientStd. errorz-statisticProb.VariableCoefficientStd. errorz-statisticProb.C0.0059190.0116520.5079940.6115C0.0063740.0102060.6245260.5323Variance equationVariance equationC4.012945739.4960.0006990.9994C0.0115820.0459430.2520910.801RESID (−1)2492.9499705091.90.0006990.9994RESID (−1)20.6417692.5009380.2566110.7975GARCH (−1)0.7586570.0865418.7664790GARCH (−1)0.7877820.115036.8484780T-DIST. DOF2.0002990.4277984.6758040T-DIST. DOF2.1486450.6431483.3408230.0008ARCH + GARCH493.708557ARCH + GARCH1.429551R-squared−0.000496Mean dependent var0.013354R-squared−0.001009Mean dependent var0.012429AdjustedR-squared−0.000496S.D. dependent var0.335083AdjustedR-squared−0.001009S.D. dependent var0.191365S.E. of regression0.335166Akaike info criterion0.157424S.E. of regression0.191462Akaike info criterion−0.736736Sum squared resid13.81739Schwarz criterion0.271145Sum squared resid4.508874Schwarz criterion−0.623015Log likelihood−4.760306Hannan–Quinn criterion.0.20362Log likelihood50.67763Hannan–Quinn criterion.−0.69054Durbin–Watson stat2.542814Durbin–Watson stat2.697083Source: author’s own work. Notes: method: ML ARCH–student’st distribution (BFGS/Marquardt steps). Sample: 7/01/2021–12/21/2021. Included observations: 124. For the definition of variables, please see Table 2.We estimated the conditional volatilities (CV) after assessing the GARCH model (1, 1), as shown in Figure7. Among the Romanian companies operating in the energy field, we notice that the highest CV was registered by the major performers in the energy allocation and supply market in Romania, respectively Electrica.Figure 7
Conditional volatility (CV) plot of Romanian stock market index and energy companies. Source: author’s own work. Notes: variables’ descriptions are provided in Table2.


#### 4.2.3. ARDL Estimation Results

Since the stationarity test confirmed that the selected variables are integrated in the order of 0 and 1, in line with Le et al. [24], this fact allows us to consider the approach of the ARDL analysis technique. This method permits the cointegration examination of variables that are stationary and non-stationary. For the proper choice of the ARDL model that would allow us to research the relationships that are established between the variables, it is imperative to choose the correct number of offsets. Therefore, alike Nyga-Łukaszewska and Aruga [39], we will analyze the Akaike Information Criteria (AIC) to select the optimal offsets for the variables included in the ARDL model. We will apply the criteria graph, which will indicate the right lags for the ARDL model, and the lowest value is preferred. Figure 8 shows the plots of criteria graph for each ARDL model which considers every energy variable (e.g., BET-NG, SNG, EL, TGN, TEL, EEX-B, EEX-P, APX, and LNGI), as well as the number of new cases at European and global level due to COVID-19. The horizontal axis of each chart represents the ARDL models estimated, and the vertical axis shows the AIC value of the models. The top 20 results are presented in the criteria graph.Figure 8
Akaike information criteria (AIC). Source: author’s own work. Notes: for the definition of variables, please see Table2.
(a)(b)The figures in Table8 signify the results for the ARDL bound test for cointegration. There are provided two critical values for the cointegration test: the lower critical bound assumes all the variables are I (0), meaning that there is no cointegration and the upper bound assumes that all the variables are I (1), meaning that there is cointegration among the variables. If the value of the F-statistic is more than the critical value of bounds, it indicates the long-run relationship between variables. Accordingly, in all cases, the value of the F-statistical test is higher than the limit of 1%, which suggests that there is a long-term relationship between variables. Hence, the null hypothesis is rejected, which means that the variables in all estimated models are cointegrated. As such, we notice cointegration relationships between COVID-19 and BET-NG, SNG, EL, TGN, TEL, EEX-B, EEX-P, APX, and LNGI.Table 8
The results of the ARDL bounds test for the model environment and COVID-19.
Test statistic:F-statisticBET-NG40.32473SNG33.30209EL39.27895TGN52.79222TEL40.24938EEX-B25.06441EEX-P25.75187APX23.92644LNGI38.06681Critical value boundsSignificanceI0 boundI1 bound10%3.174.145%3.794.852.50%4.415.521%5.156.36Source: author’s own work. Notes: null hypothesis: no long-run relationships exist. For the definition of variables, please see Table2.Table9 shows the results of the long-term relationship between variables for COVID-19 new cases in Europe and worldwide and energy market-specific variables. The quantitative outcomes provide support for no impact of the number of new cases of COVID-19 in Europe and globally on the price of energy, natural gas, the BET-NG index, and most of the companies operating in the energy field in Romania. However, there is noticed a positive long-term impact of the number of new cases in the EU and globally on TGN.Table 9
ARDL long-term term coefficients.
VariablesCoefficientStd. errort-statisticProb.APXCNE000.8006750.425D (CNW)00−0.226990.8208BET-NGCNE00−0.997160.3207D (CNW)000.3142190.7539EEX-BCNE000.6295790.5303D (CNW)00−0.438450.6619EEX-PCNE000.6713890.5034D (CNW)00−1.317880.1903ELCNE00−0.825070.411D (CNW)000.8276370.4095LNGICNE000.1109850.9118D (CNW)000.0123880.9901SNGCNE00−0.734280.4642D (CNW)001.5042020.1352TELCNE00−1.351050.1794D (CNW)000.2212260.8253TGNCNE00−2.464090.0152D (CNW)00−2.273020.0249Source: author’s own work. Notes: for the definition of variables, please see Table2.Table10 shows the short-run impact of pandemic on the energy market. For all energy market variables is registered a coefficient of the error correction term (CointEq(−1)) which is negative and significant at the 5% level of significance. Consequently, the negative and significant error correction term, which indicates the speed of conversion, exhibits that on the next day, the dependent variable will reach equilibrium with a speed of between 87% and 229%. Also, the short-term results provided in Table 10 show a positive impact of the new cases of COVID-19 registered in Europe on the evolution of the TEL share price. Therefore, an increase in the number of new cases of COVID-19 in Europe during the period under review leads to an increase in the price of TEL. Such a relationship was also identified in the case of TGN, where the new number of COVID-19 cases both in Europe and globally exerts a positive impact on the share price. Another outstanding result identified from these ARDL models is that, in the short term, the outcomes show a negative impact of new cases of COVID-19 disease in Europe on the variable EEX-B–proxy variable for the price of electricity.Table 10
ARDL cointegrating and short-term coefficients.
VariablesCoefficientStd. errort-statisticProb.Dependent variable: APXD(APX(−1))0.8090480.2277743.5519770.0006D(APX(−2))0.4349990.1630332.6681640.0087D(APX(−3))0.2530360.0916282.7615630.0067D(CNE)000.7959410.4277D(D(CNW))00−0.227210.8207CointEq(−1)−2.2976430.272947−8.417920Dependent variable: BET-NGD(CNE)00−0.993770.3224D(D(CNW))000.3146070.7536CointEq(−1)−1.0037230.09165−10.95170Dependent variable: EEX-BD(EEX-B(−1))0.8332630.1900634.3841530D(EEX-B(−2))0.4648220.1403763.3112570.0013D(EEX-B(−3))0.3138670.0876143.5823950.0005D(CNE)−0.0000010.000001−2.191950.0305D(D(CNW))00−0.439390.6612CointEq(−1)−2.0338930.235973−8.619180Dependent variable: EEX-PD(EEX-P(−1))0.8392160.2083584.0277710.0001D(EEX-P(−2))0.5149260.1524793.3770250.001D(EEX-P(−3))0.2798960.0895423.1258750.0023D(CNE)−0.0000020.000001−1.750910.0827D(D(CNW))00−0.757570.4503CointEq(−1)−2.2064870.257442−8.570810Dependent variable: ELD(CNE)00−0.824380.4114D(D(CNW))000.8330040.4065CointEq(−1)−0.990390.091465−10.82810Dependent variable: LNGID(CNE)000.1109690.9118D(D(CNW))000.0123880.9901CointEq(−1)−0.9831640.091917−10.69620Dependent variable: SNGD(CNE)00−0.737230.4625D(D(CNW))000.8354340.4052CointEq(−1)−0.8734530.089253−9.786280Dependent variable: TELD(CNE)001.3832630.1693D(CNE(−1))000.1019940.9189D(CNE(−2))002.335510.0213D(D(CNW))000.2212640.8253CointEq(−1)−1.0057320.091574−10.98270Dependent variable: TGND(CNE)002.8408560.0053D(CNE(−1))003.2864660.0013D(D(CNW))00−2.363280.0198CointEq(−1)−1.0862950.088083−12.33260Source: author’s own work. Notes: for the definition of variables, please see Table2.Further, alike Jeris and Nath [37], Geyikçi [41], this study applied cumulative sum (CUSUM) and cumulative sum of the squares (CUSUM of Squares) to confirm the stability of the long-run and short-run parameters, respectively the reliability and stability of the examined models. While the cumulative sum of squares test detects abrupt changes from the constancy of the regression coefficient, the cumulative sum test captures systematic variations in the regression coefficients [74]. Under the null hypothesis, the regression coefficients remain constant over time, being equal (or stable) in all sequential subsamples Ploberger and Krämer [75]. Hence, the estimated models are stable and there is no structural break in the observed time series, if the null hypothesis is valid [76]. Contrarily, the alternative assumption is that the regression coefficients fluctuate during the course of the sample. Nevertheless, the test is commonly questioned on the grounds that it is quite straightforward and does not request a priori information regarding the timing of the structural change [78]. Concretely, these tests produce a diagram by recurrently computing the regression coefficients and residuals [76]. For the stability of the estimated outcomes, the graph should maintain the critical values.The results of CUSUM and CUSUM of squares are plotted in Figure9. As far as the CUSUM test is concerned, all the models are stable as the CUSUM blue line is within the 5% significance boundaries illustrated by red dotted lines. Howbeit, the CUSUM test is less powerful because the confidence interval of the test is approximated [77]. In addition, the CUSUM of Squares test shows the stability of APX, BET-NG, EEX-P, and TGN models, but EEX-B, EL, LNGI, SNG, and TEL models are unstable as the CUSUM of Squares blue line crosses the 5% boundary. However, although it is known that the CUSUM of Squares test statistic is distributed as a beta random variable, there should be acknowledged that the confidence interval of the test is also approximated [77].Figure 9
CUSUM test and CUSUM of square test of each ARDL model. Source: author’s own work. Notes: for the definition of variables, please see Table2.
(a)(b)(c)Further, even though several models are unstable as pointed out by Figure9, Table 11 provides supplementary diagnostic tests in order to check for robustness. To avoid the serial correlation problem (the association among included variables and its lagged value), the Breusch-Godfrey LM test for autocorrelation (presence of autocorrelation in the null hypothesis) was conducted, while heteroscedasticity was examined by Breusch-Pagan-Godfrey test. Accordingly, the hypothesis of no serial correlation between variables and its lagged value is rejected at a 5% level of significance (Prob. Chi-Square > 0.05), except for the EEX-P model. In addition, the probability associated with the Chi-Square value is above the 0.05 significance level signifying that the assumption of homoscedasticity fails to be rejected, except EEX-B, EEX-P, and SNG models. Therefore, the errors of the EEX-B, EEX-P, and SNG models are not homoscedastic.Table 11
ARDL diagnosis tests.
APXBreusch–Godfrey serial correlation LM testF-statistic1.956852Prob.F(2, 111)0.1461Obs∗R-squared4.086932Prob. Chi-square(2)0.1296Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic0.871964Prob.F(6, 113)0.5179Obs∗R-squared5.310026Prob. Chi-square(6)0.5047Scaled explained SS13.70353Prob. Chi-square(6)0.0331BET-NGBreusch–Godfrey serial correlation LM testF-statistic0.095708Prob.F(2, 117)0.9088Obs∗R-squared0.200903Prob. Chi-square(2)0.9044Heteroskedasticity test: Breusch-Pagan-GodfreyF-statistic2.362356Prob.F(3, 119)0.0747Obs∗R-squared6.91355Prob. Chi-square(3)0.0747Scaled explained SS10.20941Prob. Chi-square(3)0.0169EEX-BBreusch–Godfrey serial correlation LM testF-statistic1.43912Prob.F(2, 110)0.2416Obs∗R-squared3.059835Prob. Chi-square(2)0.2166Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic3.030858Prob.F(7, 112)0.0059Obs∗R-squared19.11122Prob. Chi-square(7)0.0078Scaled explained SS36.21409Prob. Chi-square(7)0EEX-PBreusch–Godfrey serial correlation LM testF-statistic3.67646Prob.F(2, 109)0.0285Obs∗R-squared7.583398Prob. Chi-square(2)0.0226Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic2.361501Prob.F(8, 111)0.0219Obs∗R-squared17.45327Prob. Chi-square(8)0.0257Scaled explained SS36.76958Prob. Chi-square(8)0ELBreusch–Godfrey serial correlation LM testF-statistic0.078257Prob.F(2, 117)0.9248Obs∗R-squared0.164321Prob. Chi-square(2)0.9211Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic1.515814Prob.F(3, 119)0.214Obs∗R-squared4.527293Prob. Chi-square(3)0.2099Scaled explained SS8.60173Prob. Chi-square(3)0.0351LNGIBreusch–Godfrey serial correlation LM testF-statistic0.03895Prob.F(2, 117)0.9618Obs∗R-squared0.081841Prob. Chi-square(2)0.9599Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic0.019132Prob.F(3, 119)0.9964Obs∗R-squared0.059298Prob. Chi-square(3)0.9962Scaled explained SS1.165012Prob. Chi-square(3)0.7614SNGBreusch–Godfrey serial correlation LM testF-statistic0.210536Prob.F(2, 115)0.8105Obs∗R-squared0.445074Prob. Chi-square(2)0.8005Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic4.599008Prob.F(4, 117)0.0017Obs∗R-squared16.57594Prob. Chi-square(4)0.0023Scaled explained SS38.89259Prob. Chi-square(4)0TELBreusch–Godfrey serial correlation LM testF-statistic0.4525Prob.F(2, 112)0.6372Obs∗R-squared0.969886Prob. Chi-square(2)0.6157Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic0.837925Prob.F(6, 114)0.5431Obs∗R-squared5.110863Prob. Chi-square(6)0.5297Scaled explained SS14.91187Prob. Chi-square(6)0.021TGNBreusch–Godfrey serial correlation LM testF-statistic0.246637Prob.F(2, 114)0.7818Obs∗R-squared0.525616Prob. Chi-square(2)0.7689Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic1.580195Prob.F(5, 116)0.1711Obs∗R-squared7.779751Prob. Chi-square(5)0.1688Scaled explained SS10.33263Prob. Chi-square(5)0.0663Source: author’s own work. Notes: for the definition of variables, please see Table2.


#### 4.2.4. Variance Decomposition Research

To determine the extent to which COVID-19 variables contribute to the explanation of energy market variables, the variance decomposition approach was used, the outcomes being revealed in Table12. Decomposition variation indicates the extent to which a certain variable can explain the evolution of the variation of another variable. In addition, it shows which of the independent variables is stronger in explaining the variability of dependent variables over time.Table 12
The results of the variance decomposition.
Variance decomposition of BET-NGPeriodBET-NGCNED(CNW)498.297860.5159061.186235Variance decomposition of APXPeriodAPXCNED(CNW)495.799620.5795113.620871Variance decomposition of EEX-BPeriodEEX-BCNED(CNW)498.267531.3516390.380828Variance decomposition of EEX-PPeriodEEX-PCNED(CNW)497.445531.8299810.724489Variance decomposition of LNGIPeriodLNGICNED(CNW)497.07622.1712360.75256Variance decomposition of ELPeriodELCNED(CNW)497.454331.6309740.914701Variance decomposition of SNGPeriodSNGCNED(CNW)496.267662.227231.505114Variance decomposition of TELPeriodTELCNED(CNW)498.295610.7208730.983517Variance decomposition of TGNPeriodTGNCNED(CNW)491.53255.3889153.078584Source: author’s own work. Notes: null hypothesis: no long-run relationships exist. Notes: for the definition of variables, please see Table2.From the results obtained after the decomposition of the variance, it can be noticed that the number of new cases of COVID-19 registered in Europe shows a higher contribution in the case of EEX-B, EEX-P, LNGI, EL, SNG, and TGN, whereas for the variables BET-NG, APX, TEL, the largest contribution comes from the number of new COVID-19 cases registered globally.


#### 4.2.5. Causality Assessment

To explore causality between selected variables, the Granger causality test is applied. To employ the Granger causality test, the data series must be stationary and have therefore been converted into stationary series. Table13 shows the results after the Granger causality test for the energy market and COVID-19 variables.Table 13
Results of granger causality test.
Null hypothesisF-statisticProb.D(D(CNW)) does not granger cause D(CNE)0.715270.6381D(CNE) does not granger cause D(D(CNW))1.34730.2432APX does not granger cause D(CNE)1.0480.3989D(CNE) does not granger cause APX0.556420.764BET-NG does not granger cause D(CNE)1.013430.4207D(CNE) does not granger cause BET-NG1.564390.1649EEX-P does not granger cause D(CNE)1.561610.1658D(CNE) does not granger cause EEX-P1.454460.2012EEX-B does not granger cause D(CNE)1.856280.0954D(CNE) does not granger cause EEX-B1.261170.2818EL does not granger cause D(CNE)1.020950.4159D(CNE) does not granger cause EL0.858160.5284SNG does not granger cause D(CNE)1.021230.4157D(CNE) does not granger cause SNG1.194950.3149TEL does not granger cause D(CNE)0.868930.5204D(CNE) does not granger cause TEL1.992310.0733TGN does not granger cause D(CNE)1.213910.3051D(CNE) does not granger cause TGN3.895490.0015APX does not granger cause D(D(CNW))0.353650.9063D(D(CNW)) does not granger cause APX0.129340.9924BET-NG does not granger cause D(D(CNW))1.948860.0799D(D(CNW)) does not granger cause BET-NG1.761130.1144EEX-P does not granger cause D(D(CNW))0.435350.8537D(D(CNW)) does not granger cause EEX-P0.663650.6791EEX-B does not granger cause D(D(CNW))0.32820.9208D(D(CNW)) does not granger cause EEX-B0.487020.8167EL does not granger cause D(D(CNW))0.584540.7419D(D(CNW)) does not granger cause EL0.43870.8514SNG does not granger cause D(D(CNW))1.84630.0973D(D(CNW)) does not granger cause SNG1.719580.1238TEL does not granger cause D(D(CNW))1.932760.0824D(D(CNW)) does not granger cause TEL0.969730.4496TGN does not granger cause D(D(CNW))0.767930.5968D(D(CNW)) does not granger cause TGN1.760570.1146Source: author’s own work. Notes: sample: 7/01/2021–12/21/2021. For the definition of variables, please see Table2.The causality test confirms that between the number of new cases of COVID-19 registered in Europe and TEL, respectively TGN, there is a one-way relationship running from CNE to both energy companies listed on the BSE. From an econometric point of view, we came to this decision because thep-value is below the threshold of 10% and 5% respectively, which leads to the rejection of the null hypothesis: the number of new cases of COVID-19 in Europe does not determine Granger-type causation variables analyzed.To our knowledge, the existing studies on the impact of COVID-19 on the energy market during July-December 2021, did not address this type of relationship study, the novelty of the article deepening consisting in this ARDL model that will allow the analysis of long-term relationships between variables selected, as well as the decomposition of the variance and the identification of the causal relations.The innovations of this research include the following aspects. First, this study measured the impact of the COVID-19 pandemic on the energy market. Second, this research integrated advanced econometric models to obtain detailed results on the long-term or short-term relationships between COVID-19 variables and the energy market. Third, this study will help fill the gap in the literature and will be a focal point for future energy market research during the COVID-19 pandemic.


## 4.1. Summary Statistics

The descriptive statistics of the variables are provided in Table3. The skewness and kurtosis indicators indicate the deviation in relation to a symmetric distribution around the average, thus suggesting the degree of flattening or sharpening. A kurtosis greater than the value of three implies that the returns of the indices show heavy tails than the normal distribution. Specifically, the probability of extreme returns is higher than the probability that they are below normal distribution. This feature is called leptokurtic or basically heavy tails. A positive skewness signifies an asymmetric distribution on the right and a negative skewness on the left. For a series with a normal distribution, kurtosis takes the value of three. However, for a value less than three the distribution is flatter than the normal one (e.g., platykurtic), whereas for kurtosis greater than three, the distribution is leptokurtic. According to Table 3, most of the variables exhibit a value of the kurtosis greater than three, thus presenting a leptokurtic distribution, consistent with Bildirici et al. [8], Wu and Ma [9], Wang et al. [36], Geyikçi [41], Bourghelle et al. [48], Dutta et al. [53]. Consequently, extreme negative returns are much more likely to occur than normal distribution forecasts.Table 3
Descriptive statistics of the variables.
VariablesMeanStd. dev.SkewnessKurtosisJarque-BeraProbabilityCNE2310541089950.892.4418.010.00CNW5788641130850.352.384.500.11BET-NG0.00010.01−0.744.4021.550.00SNG0.00120.01−0.826.5679.630.00EL−0.00200.01−0.755.2337.410.00TEL−0.00120.020.488.58165.980.00TGN−0.00220.01−0.374.5214.750.00EEX-B0.01240.270.155.4932.580.00EEX-P0.01340.340.105.8843.040.00APX0.01240.190.596.5371.430.00LNGI0.00990.062.0643.218442.140.00Source: author’s own work. Notes: for the definition of variables, please see Table2.The normality of the variable distribution is provided by the Jarque-Bera test. Table3 also shows the results of the Jarque-Bera test, which indicate that the distribution of the variables is not distributed normally, in line with Villarreal-Samaniego [2], Bildirici et al. [8], Ali et al. [25], Si et al. [28], Wang et al. [36], Mensi et al. [38], Akhtaruzzaman et al. [69], Gharib et al. [71], Khan et al. [72]. Because the Jarque-Bera statistic is significant (except for the number of new cases of COVID-19 registered globally), we reject the null hypothesis of normality. The test values are quite far from the corresponding normal distribution, which supports that the series is not normally distributed. Figure 1 exhibits the density, distribution, and quantile-quantile (QQ) plots. In line with the outcomes reported in Table 3, the distribution of the selected series is dissimilar from the normal one.Figure 1
Density, distribution, and QQ plots for daily returns. Source: author’s own work. Notes: variables’ descriptions are provided in Table2.The correlations between included variables are plotted in Figure2. Alike [37], the correlations were reduced. Therefore, our empirical outcomes will not be affected by the multicollinearity issue.Figure 2
Correlations among selected variables. Source: author’s own work. Notes: variables’ descriptions are provided in Table2.Figure3 shows the evolution of the number of new cases at the European and global levels due to COVID-19. Cássaro and Pires [73] argued that the number of cases is growing rapidly, which has been achieved by stability later, this mode is called a step function.Figure 3
Evolution of COVID-19 pandemic variables. Source: author’s own work. Notes: for the definition of variables, please see Table2.Further, the daily evolution of the BET-NG index and the companies listed on BSE are presented in Figure4. During the third and fourth quarters of 2021, BET-NG, the index of the 10 Romanian companies in the energy and utility sectors, registered a relatively steady evolution. EL, SNG, and TGN, which are among the most traded companies on BSE, had a similar evolution to the BET-NG index, but TEL was marked by slightly more significant episodes of volatility than the rest of the energy companies.Figure 4
Evolution of BET-NG index and of the companies listed on BSE. Source: author’s own work. Notes: for the definition of variables, please see Table2.Likewise, the evolution related to the return of the Physical Electricity Index (Phelix), Amsterdam Power Exchange (APX), and London Natural Gas Index is exhibited in Figure5. Energy prices have witnessed high episodes of volatility compared to the price of natural gas. The return of global activity, as the restrictions imposed by the pandemic and the complete reopening of economies have lifted, has led to a rapid increase in demand for natural gas, both for electricity production and for manufacturing reasons.Figure 5
Evolution of energy market. Source: author’s own work. Notes: for the definition of variables, please see Table2.


## 4.2. The Outcomes of Time Series Investigation

### 4.2.1. Stationarity Analysis

Tables4 and 5 show the outcomes of stationarity analysis for series at the level and in the first difference through the augmented Dickey–Fuller and Kwiatkowski–Phillips–Schmidt–Shin unit root test. A couple of variables are stationary at the first difference (e.g., new cases of COVID-19 worldwide and in Europe) because the probability is above the 1% and 5% relevance level, while some are already stationary (e.g., BET-NG, SNG, EL, TGN, TEL, EEX-B, EEX-P, APX, and LNGI). Alike Bildirici et al. [8], Ali et al. [25], Si et al. [28], Mensi et al. [38], and Yousaf [47], we notice a common integration order of I (0), exception makes the indicators related to the evolution of COVID-19, so in this case, we can reject the null hypothesis and conclude that the series is not stationary.Table 4
The outcomes of the augmented Dickey–Fuller unit root test.
VariablesLevel1st differencet-statisticProb.t-statisticProb.BET-NG−10.968150.00000−8.025040.00000SNG−9.797470.00000−7.022060.00000EL−10.892070.00000−8.212170.00000TGN−11.573740.00000−8.255440.00000TEL−10.629240.00000−7.738480.00000EEX-B−9.284200.00000−10.038590.00000EEX-P−9.262510.00000−7.880230.00000APX−3.852240.00330−9.075060.00000LNGI−10.804120.00000−6.786370.00000D (CNW)−1.054550.73150−6.464690.00000CNE1.492120.99920−3.174390.02410Source: author’s own work. Notes: null hypothesis: each series has a unit root. Intercept included in test equation. Lag length: automatic selection based on Akaike info criterion. Test critical values: 1% level: −3.484198; 5% level: −2.885051; 10% level: −2.579386. For the definition of variables, please see Table2.Table 5
The outcomes of the Kwiatkowski–Phillips–Schmidt–Shin unit root test.
VariablesLevel1st differenceLM-stat.LM-stat.BET-NG0.123170.10802SNG0.051320.27447EL0.073580.22224TGN0.335890.09240TEL0.095700.03197EEX-B0.264770.10671EEX-P0.254600.10521APX0.300140.14376LNGI0.068380.33024D (CNW)0.119720.09233CNE1.172860.13543Source: author’s own work. Notes: null hypothesis: each series is stationary. Intercept included in test equation. Asymptotic critical values: 1% level: 0.739000; 5% level: 0.463000; 10% level: 0.347000. For the definition of variables, please see Table2.Further, Table6 reveals the outcomes of stationarity analysis of the variables with one structural break by means of the Zivot-Andrews unit root test, alike Zhang et al. [7], Le et al. [24]. Accordingly, the occurrence of a structural break in our sample is proved, while the mix integration is strengthened.Table 6
The outcomes of the Zivot–Andrews unit root test.
VariablesLevel1st differencet-statisticProb.Chosen breakpointt-statisticProb.Chosen breakpointBET-NG−11.602500.0145810/29/2021−10.169530.369439/28/2021SNG−10.093280.1099810/28/2021−8.848800.3967611/11/2021EL−11.146200.0322710/28/2021−7.832430.3128411/22/2021TGN−12.346240.0098311/05/2021−9.038700.1122611/24/2021TEL−10.828270.0950610/26/2021−8.099400.0481410/22/2021EEX-B−9.512800.0407710/25/2021−9.625170.1805010/22/2021EEX-P−9.478390.2147810/11/2021−11.075400.1708110/22/2021APX−8.482670.0356010/14/2021−10.477770.0891710/08/2021LNGI−11.552590.0070910/07/2021−8.209490.0262510/04/2021D (CNW)−9.573620.000028/23/2021−15.985630.342039/16/2021CNE−2.455740.0013010/19/2021−6.629530.0239310/06/2021Source: author’s own work. Notes: null hypothesis: each series has a unit root with a structural break in the intercept. Break included in the intercept. 1% critical value: −5.34.5% critical value: −4.93.10% critical value: −4.58. For the definition of variables, please see Table2.


### 4.2.2. Volatility Examination

Because all Jarque–Bera figures are larger than three as Bildirici et al. [8], Zhang and Wang [57] found, and the distributions of the selected series diverge from normally distributed data; GARCH models are required and appropriate for highlighting the progress of volatility. Further, the autocorrelation function (ACF) and partial autocorrelation function (PACF) are used to test for heteroskedasticity. In this regard, Figure 6 exhibits the plots of ACF and PACF. Thus, the existence of serial correlation is inherent in the correlogram of squared returns, except LNGI for which the GARCH model will not be estimated.Figure 6
Plots of autocorrelation function (ACF) and partial autocorrelation function (PACF) for squared returns. Source: author’s own work. Notes: the lag length was set to 20. Variables’ descriptions are provided in Table2.Table7 reports the findings after estimating the GARCH (1, 1) model. In the case of TEL, the sum of ARCH and GARCH parameters is roughly equal to one (e.g., 0.989732), in line with Napon and Asama [4], Bildirici et al. [8]. Hence, at time t, a shock wave will last longer. Specifically, the conditional variance is persistent consistent with Maneejuk et al. [44], Iglesias, and Rivera-Alonso [51]. The findings, however, show the existence of a mean-reverting process because the sum of the ARCH and GARCH effects is less than one. However, in the case of EEX-B, EEX-P, and APX, the sum of parameters is larger than one, thus suggesting that the conditional variance process is explosive.Table 7
The outcomes of the GARCH (1, 1) model.
Dependent variable: BET-NGDependent variable: SNGVariableCoefficientStd. errorz-statisticProb.VariableCoefficientStd. errorz-statisticProb.C0.0005120.0007770.6589560.5099C0.0022260.0009642.307890.021Variance equationVariance equationC4.21E − 054.80E − 050.8768020.3806C9.67E − 050.0001810.5333310.5938RESID (−1)2−0.074930.068867−1.0880460.2766RESID (−1)2−0.048840.072398−0.6746020.4999GARCH (−1)0.6563340.4563521.438220.1504GARCH (−1)0.5420730.9169320.5911810.5544T-DIST. DOF5.0615213.1427831.6105220.1073T-DIST. DOF3.4526231.2507472.7604480.0058ARCH + GARCH0.581404ARCH + GARCH0.493233R-squared−0.001889Mean dependent var8.42E − 05R-squared−0.006586Mean dependent var0.001153AdjustedR-squared−0.001889S.D. dependent var0.009884AdjustedR-squared−0.006586S.D. dependent var0.013269S.E. of regression0.009894Akaike info criterion−6.409628S.E. of regression0.013313Akaike info criterion−5.904976Sum squared resid0.01204Schwarz criterion−6.295907Sum squared resid0.0218Schwarz criterion−5.791255Log likelihood402.3969Hannan–Quinn criterion.−6.363432Log likelihood371.1085Hannan–Quinn criterion.−5.85878Durbin-Watson stat1.988924Durbin-Watson stat1.7568Dependent variable: ELDependent variable: TGNVariableCoefficientStd. errorz-statisticProb.VariableCoefficientStd. errorz-statisticProb.C−8.80E − 060.000686−0.0128180.9898C−0.0015810.000705−2.2426210.0249Variance equationVariance equationC0.0003490.0082760.0421860.9664C6.14E − 065.86E − 061.0488260.2943RESID (−1)2−1.78439542.03255−0.0424530.9661RESID (−1)20.0757880.0746651.0150350.3101GARCH (−1)0.9877170.01889952.262890GARCH (−1)0.8743390.1022928.5474640T-DIST. DOF2.0243890.584223.4651160.0005T-DIST. DOF4.3495662.4443851.7794120.0752ARCH + GARCH−0.796678ARCH + GARCH0.950127R-squared−0.025806Mean dependent var−0.001991R-squared−0.004356Mean dependent var−0.002206AdjustedR-squared−0.025806S.D. dependent var0.012388AdjustedR-squared−0.004356S.D. dependent var0.009498S.E. of regression0.012547Akaike info criterion−6.171115S.E. of regression0.009519Akaike info criterion−6.559459Sum squared resid0.019363Schwarz criterion−6.057394Sum squared resid0.011145Schwarz criterion−6.445738Log likelihood387.6091Hannan–Quinn criterion.−6.124919Log likelihood411.6864Hannan–Quinn criterion.−6.513262Durbin–Watson stat1.93045Durbin–Watson stat2.091561Dependent variable: TELDependent variable: EEX-BVariableCoefficientStd. errorz-statisticProb.VariableCoefficientStd. errorz-statisticProb.C−0.0008560.000926−0.9240740.3554C0.0061580.0107720.5716660.5675Variance equationVariance equationC5.27E − 061.80E − 062.9369290.0033C0.0060220.0147580.4080330.6832RESID (−1)2−0.050780.016305−3.1144110.0018RESID (−1)21.519244.012510.3786260.705GARCH (−1)1.0405120.0139374.696550GARCH (−1)0.671250.1101996.0912320T-DIST. DOF5.7675783.2191841.7916270.0732T-DIST. DOF2.2005450.5998483.6685030.0002ARCH + GARCH0.989732ARCH + GARCH2.19049R-squared−0.000593Mean dependent var−0.001222R-squared−0.000556Mean dependent var0.012439AdjustedR-squared−0.000593S.D. dependent var0.015068AdjustedR-squared−0.000556S.D. dependent var0.267469S.E. of regression0.015073Akaike info criterion−5.758566S.E. of regression0.267543Akaike info criterion−0.16669Sum squared resid0.027944Schwarz criterion−5.644845Sum squared resid8.804247Schwarz criterion−0.052969Log likelihood362.0311Hannan–Quinn criterion.−5.71237Log likelihood15.3348Hannan–Quinn criterion.−0.120494Durbin-Watson stat1.930154Durbin-Watson stat2.342452Dependent variable: EEX-PDependent variable: APXVariableCoefficientStd. errorz-statisticProb.VariableCoefficientStd. errorz-statisticProb.C0.0059190.0116520.5079940.6115C0.0063740.0102060.6245260.5323Variance equationVariance equationC4.012945739.4960.0006990.9994C0.0115820.0459430.2520910.801RESID (−1)2492.9499705091.90.0006990.9994RESID (−1)20.6417692.5009380.2566110.7975GARCH (−1)0.7586570.0865418.7664790GARCH (−1)0.7877820.115036.8484780T-DIST. DOF2.0002990.4277984.6758040T-DIST. DOF2.1486450.6431483.3408230.0008ARCH + GARCH493.708557ARCH + GARCH1.429551R-squared−0.000496Mean dependent var0.013354R-squared−0.001009Mean dependent var0.012429AdjustedR-squared−0.000496S.D. dependent var0.335083AdjustedR-squared−0.001009S.D. dependent var0.191365S.E. of regression0.335166Akaike info criterion0.157424S.E. of regression0.191462Akaike info criterion−0.736736Sum squared resid13.81739Schwarz criterion0.271145Sum squared resid4.508874Schwarz criterion−0.623015Log likelihood−4.760306Hannan–Quinn criterion.0.20362Log likelihood50.67763Hannan–Quinn criterion.−0.69054Durbin–Watson stat2.542814Durbin–Watson stat2.697083Source: author’s own work. Notes: method: ML ARCH–student’st distribution (BFGS/Marquardt steps). Sample: 7/01/2021–12/21/2021. Included observations: 124. For the definition of variables, please see Table 2.We estimated the conditional volatilities (CV) after assessing the GARCH model (1, 1), as shown in Figure7. Among the Romanian companies operating in the energy field, we notice that the highest CV was registered by the major performers in the energy allocation and supply market in Romania, respectively Electrica.Figure 7
Conditional volatility (CV) plot of Romanian stock market index and energy companies. Source: author’s own work. Notes: variables’ descriptions are provided in Table2.


### 4.2.3. ARDL Estimation Results

Since the stationarity test confirmed that the selected variables are integrated in the order of 0 and 1, in line with Le et al. [24], this fact allows us to consider the approach of the ARDL analysis technique. This method permits the cointegration examination of variables that are stationary and non-stationary. For the proper choice of the ARDL model that would allow us to research the relationships that are established between the variables, it is imperative to choose the correct number of offsets. Therefore, alike Nyga-Łukaszewska and Aruga [39], we will analyze the Akaike Information Criteria (AIC) to select the optimal offsets for the variables included in the ARDL model. We will apply the criteria graph, which will indicate the right lags for the ARDL model, and the lowest value is preferred. Figure 8 shows the plots of criteria graph for each ARDL model which considers every energy variable (e.g., BET-NG, SNG, EL, TGN, TEL, EEX-B, EEX-P, APX, and LNGI), as well as the number of new cases at European and global level due to COVID-19. The horizontal axis of each chart represents the ARDL models estimated, and the vertical axis shows the AIC value of the models. The top 20 results are presented in the criteria graph.Figure 8
Akaike information criteria (AIC). Source: author’s own work. Notes: for the definition of variables, please see Table2.
(a)(b)The figures in Table8 signify the results for the ARDL bound test for cointegration. There are provided two critical values for the cointegration test: the lower critical bound assumes all the variables are I (0), meaning that there is no cointegration and the upper bound assumes that all the variables are I (1), meaning that there is cointegration among the variables. If the value of the F-statistic is more than the critical value of bounds, it indicates the long-run relationship between variables. Accordingly, in all cases, the value of the F-statistical test is higher than the limit of 1%, which suggests that there is a long-term relationship between variables. Hence, the null hypothesis is rejected, which means that the variables in all estimated models are cointegrated. As such, we notice cointegration relationships between COVID-19 and BET-NG, SNG, EL, TGN, TEL, EEX-B, EEX-P, APX, and LNGI.Table 8
The results of the ARDL bounds test for the model environment and COVID-19.
Test statistic:F-statisticBET-NG40.32473SNG33.30209EL39.27895TGN52.79222TEL40.24938EEX-B25.06441EEX-P25.75187APX23.92644LNGI38.06681Critical value boundsSignificanceI0 boundI1 bound10%3.174.145%3.794.852.50%4.415.521%5.156.36Source: author’s own work. Notes: null hypothesis: no long-run relationships exist. For the definition of variables, please see Table2.Table9 shows the results of the long-term relationship between variables for COVID-19 new cases in Europe and worldwide and energy market-specific variables. The quantitative outcomes provide support for no impact of the number of new cases of COVID-19 in Europe and globally on the price of energy, natural gas, the BET-NG index, and most of the companies operating in the energy field in Romania. However, there is noticed a positive long-term impact of the number of new cases in the EU and globally on TGN.Table 9
ARDL long-term term coefficients.
VariablesCoefficientStd. errort-statisticProb.APXCNE000.8006750.425D (CNW)00−0.226990.8208BET-NGCNE00−0.997160.3207D (CNW)000.3142190.7539EEX-BCNE000.6295790.5303D (CNW)00−0.438450.6619EEX-PCNE000.6713890.5034D (CNW)00−1.317880.1903ELCNE00−0.825070.411D (CNW)000.8276370.4095LNGICNE000.1109850.9118D (CNW)000.0123880.9901SNGCNE00−0.734280.4642D (CNW)001.5042020.1352TELCNE00−1.351050.1794D (CNW)000.2212260.8253TGNCNE00−2.464090.0152D (CNW)00−2.273020.0249Source: author’s own work. Notes: for the definition of variables, please see Table2.Table10 shows the short-run impact of pandemic on the energy market. For all energy market variables is registered a coefficient of the error correction term (CointEq(−1)) which is negative and significant at the 5% level of significance. Consequently, the negative and significant error correction term, which indicates the speed of conversion, exhibits that on the next day, the dependent variable will reach equilibrium with a speed of between 87% and 229%. Also, the short-term results provided in Table 10 show a positive impact of the new cases of COVID-19 registered in Europe on the evolution of the TEL share price. Therefore, an increase in the number of new cases of COVID-19 in Europe during the period under review leads to an increase in the price of TEL. Such a relationship was also identified in the case of TGN, where the new number of COVID-19 cases both in Europe and globally exerts a positive impact on the share price. Another outstanding result identified from these ARDL models is that, in the short term, the outcomes show a negative impact of new cases of COVID-19 disease in Europe on the variable EEX-B–proxy variable for the price of electricity.Table 10
ARDL cointegrating and short-term coefficients.
VariablesCoefficientStd. errort-statisticProb.Dependent variable: APXD(APX(−1))0.8090480.2277743.5519770.0006D(APX(−2))0.4349990.1630332.6681640.0087D(APX(−3))0.2530360.0916282.7615630.0067D(CNE)000.7959410.4277D(D(CNW))00−0.227210.8207CointEq(−1)−2.2976430.272947−8.417920Dependent variable: BET-NGD(CNE)00−0.993770.3224D(D(CNW))000.3146070.7536CointEq(−1)−1.0037230.09165−10.95170Dependent variable: EEX-BD(EEX-B(−1))0.8332630.1900634.3841530D(EEX-B(−2))0.4648220.1403763.3112570.0013D(EEX-B(−3))0.3138670.0876143.5823950.0005D(CNE)−0.0000010.000001−2.191950.0305D(D(CNW))00−0.439390.6612CointEq(−1)−2.0338930.235973−8.619180Dependent variable: EEX-PD(EEX-P(−1))0.8392160.2083584.0277710.0001D(EEX-P(−2))0.5149260.1524793.3770250.001D(EEX-P(−3))0.2798960.0895423.1258750.0023D(CNE)−0.0000020.000001−1.750910.0827D(D(CNW))00−0.757570.4503CointEq(−1)−2.2064870.257442−8.570810Dependent variable: ELD(CNE)00−0.824380.4114D(D(CNW))000.8330040.4065CointEq(−1)−0.990390.091465−10.82810Dependent variable: LNGID(CNE)000.1109690.9118D(D(CNW))000.0123880.9901CointEq(−1)−0.9831640.091917−10.69620Dependent variable: SNGD(CNE)00−0.737230.4625D(D(CNW))000.8354340.4052CointEq(−1)−0.8734530.089253−9.786280Dependent variable: TELD(CNE)001.3832630.1693D(CNE(−1))000.1019940.9189D(CNE(−2))002.335510.0213D(D(CNW))000.2212640.8253CointEq(−1)−1.0057320.091574−10.98270Dependent variable: TGND(CNE)002.8408560.0053D(CNE(−1))003.2864660.0013D(D(CNW))00−2.363280.0198CointEq(−1)−1.0862950.088083−12.33260Source: author’s own work. Notes: for the definition of variables, please see Table2.Further, alike Jeris and Nath [37], Geyikçi [41], this study applied cumulative sum (CUSUM) and cumulative sum of the squares (CUSUM of Squares) to confirm the stability of the long-run and short-run parameters, respectively the reliability and stability of the examined models. While the cumulative sum of squares test detects abrupt changes from the constancy of the regression coefficient, the cumulative sum test captures systematic variations in the regression coefficients [74]. Under the null hypothesis, the regression coefficients remain constant over time, being equal (or stable) in all sequential subsamples Ploberger and Krämer [75]. Hence, the estimated models are stable and there is no structural break in the observed time series, if the null hypothesis is valid [76]. Contrarily, the alternative assumption is that the regression coefficients fluctuate during the course of the sample. Nevertheless, the test is commonly questioned on the grounds that it is quite straightforward and does not request a priori information regarding the timing of the structural change [78]. Concretely, these tests produce a diagram by recurrently computing the regression coefficients and residuals [76]. For the stability of the estimated outcomes, the graph should maintain the critical values.The results of CUSUM and CUSUM of squares are plotted in Figure9. As far as the CUSUM test is concerned, all the models are stable as the CUSUM blue line is within the 5% significance boundaries illustrated by red dotted lines. Howbeit, the CUSUM test is less powerful because the confidence interval of the test is approximated [77]. In addition, the CUSUM of Squares test shows the stability of APX, BET-NG, EEX-P, and TGN models, but EEX-B, EL, LNGI, SNG, and TEL models are unstable as the CUSUM of Squares blue line crosses the 5% boundary. However, although it is known that the CUSUM of Squares test statistic is distributed as a beta random variable, there should be acknowledged that the confidence interval of the test is also approximated [77].Figure 9
CUSUM test and CUSUM of square test of each ARDL model. Source: author’s own work. Notes: for the definition of variables, please see Table2.
(a)(b)(c)Further, even though several models are unstable as pointed out by Figure9, Table 11 provides supplementary diagnostic tests in order to check for robustness. To avoid the serial correlation problem (the association among included variables and its lagged value), the Breusch-Godfrey LM test for autocorrelation (presence of autocorrelation in the null hypothesis) was conducted, while heteroscedasticity was examined by Breusch-Pagan-Godfrey test. Accordingly, the hypothesis of no serial correlation between variables and its lagged value is rejected at a 5% level of significance (Prob. Chi-Square > 0.05), except for the EEX-P model. In addition, the probability associated with the Chi-Square value is above the 0.05 significance level signifying that the assumption of homoscedasticity fails to be rejected, except EEX-B, EEX-P, and SNG models. Therefore, the errors of the EEX-B, EEX-P, and SNG models are not homoscedastic.Table 11
ARDL diagnosis tests.
APXBreusch–Godfrey serial correlation LM testF-statistic1.956852Prob.F(2, 111)0.1461Obs∗R-squared4.086932Prob. Chi-square(2)0.1296Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic0.871964Prob.F(6, 113)0.5179Obs∗R-squared5.310026Prob. Chi-square(6)0.5047Scaled explained SS13.70353Prob. Chi-square(6)0.0331BET-NGBreusch–Godfrey serial correlation LM testF-statistic0.095708Prob.F(2, 117)0.9088Obs∗R-squared0.200903Prob. Chi-square(2)0.9044Heteroskedasticity test: Breusch-Pagan-GodfreyF-statistic2.362356Prob.F(3, 119)0.0747Obs∗R-squared6.91355Prob. Chi-square(3)0.0747Scaled explained SS10.20941Prob. Chi-square(3)0.0169EEX-BBreusch–Godfrey serial correlation LM testF-statistic1.43912Prob.F(2, 110)0.2416Obs∗R-squared3.059835Prob. Chi-square(2)0.2166Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic3.030858Prob.F(7, 112)0.0059Obs∗R-squared19.11122Prob. Chi-square(7)0.0078Scaled explained SS36.21409Prob. Chi-square(7)0EEX-PBreusch–Godfrey serial correlation LM testF-statistic3.67646Prob.F(2, 109)0.0285Obs∗R-squared7.583398Prob. Chi-square(2)0.0226Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic2.361501Prob.F(8, 111)0.0219Obs∗R-squared17.45327Prob. Chi-square(8)0.0257Scaled explained SS36.76958Prob. Chi-square(8)0ELBreusch–Godfrey serial correlation LM testF-statistic0.078257Prob.F(2, 117)0.9248Obs∗R-squared0.164321Prob. Chi-square(2)0.9211Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic1.515814Prob.F(3, 119)0.214Obs∗R-squared4.527293Prob. Chi-square(3)0.2099Scaled explained SS8.60173Prob. Chi-square(3)0.0351LNGIBreusch–Godfrey serial correlation LM testF-statistic0.03895Prob.F(2, 117)0.9618Obs∗R-squared0.081841Prob. Chi-square(2)0.9599Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic0.019132Prob.F(3, 119)0.9964Obs∗R-squared0.059298Prob. Chi-square(3)0.9962Scaled explained SS1.165012Prob. Chi-square(3)0.7614SNGBreusch–Godfrey serial correlation LM testF-statistic0.210536Prob.F(2, 115)0.8105Obs∗R-squared0.445074Prob. Chi-square(2)0.8005Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic4.599008Prob.F(4, 117)0.0017Obs∗R-squared16.57594Prob. Chi-square(4)0.0023Scaled explained SS38.89259Prob. Chi-square(4)0TELBreusch–Godfrey serial correlation LM testF-statistic0.4525Prob.F(2, 112)0.6372Obs∗R-squared0.969886Prob. Chi-square(2)0.6157Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic0.837925Prob.F(6, 114)0.5431Obs∗R-squared5.110863Prob. Chi-square(6)0.5297Scaled explained SS14.91187Prob. Chi-square(6)0.021TGNBreusch–Godfrey serial correlation LM testF-statistic0.246637Prob.F(2, 114)0.7818Obs∗R-squared0.525616Prob. Chi-square(2)0.7689Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic1.580195Prob.F(5, 116)0.1711Obs∗R-squared7.779751Prob. Chi-square(5)0.1688Scaled explained SS10.33263Prob. Chi-square(5)0.0663Source: author’s own work. Notes: for the definition of variables, please see Table2.


### 4.2.4. Variance Decomposition Research

To determine the extent to which COVID-19 variables contribute to the explanation of energy market variables, the variance decomposition approach was used, the outcomes being revealed in Table12. Decomposition variation indicates the extent to which a certain variable can explain the evolution of the variation of another variable. In addition, it shows which of the independent variables is stronger in explaining the variability of dependent variables over time.Table 12
The results of the variance decomposition.
Variance decomposition of BET-NGPeriodBET-NGCNED(CNW)498.297860.5159061.186235Variance decomposition of APXPeriodAPXCNED(CNW)495.799620.5795113.620871Variance decomposition of EEX-BPeriodEEX-BCNED(CNW)498.267531.3516390.380828Variance decomposition of EEX-PPeriodEEX-PCNED(CNW)497.445531.8299810.724489Variance decomposition of LNGIPeriodLNGICNED(CNW)497.07622.1712360.75256Variance decomposition of ELPeriodELCNED(CNW)497.454331.6309740.914701Variance decomposition of SNGPeriodSNGCNED(CNW)496.267662.227231.505114Variance decomposition of TELPeriodTELCNED(CNW)498.295610.7208730.983517Variance decomposition of TGNPeriodTGNCNED(CNW)491.53255.3889153.078584Source: author’s own work. Notes: null hypothesis: no long-run relationships exist. Notes: for the definition of variables, please see Table2.From the results obtained after the decomposition of the variance, it can be noticed that the number of new cases of COVID-19 registered in Europe shows a higher contribution in the case of EEX-B, EEX-P, LNGI, EL, SNG, and TGN, whereas for the variables BET-NG, APX, TEL, the largest contribution comes from the number of new COVID-19 cases registered globally.


### 4.2.5. Causality Assessment

To explore causality between selected variables, the Granger causality test is applied. To employ the Granger causality test, the data series must be stationary and have therefore been converted into stationary series. Table13 shows the results after the Granger causality test for the energy market and COVID-19 variables.Table 13
Results of granger causality test.
Null hypothesisF-statisticProb.D(D(CNW)) does not granger cause D(CNE)0.715270.6381D(CNE) does not granger cause D(D(CNW))1.34730.2432APX does not granger cause D(CNE)1.0480.3989D(CNE) does not granger cause APX0.556420.764BET-NG does not granger cause D(CNE)1.013430.4207D(CNE) does not granger cause BET-NG1.564390.1649EEX-P does not granger cause D(CNE)1.561610.1658D(CNE) does not granger cause EEX-P1.454460.2012EEX-B does not granger cause D(CNE)1.856280.0954D(CNE) does not granger cause EEX-B1.261170.2818EL does not granger cause D(CNE)1.020950.4159D(CNE) does not granger cause EL0.858160.5284SNG does not granger cause D(CNE)1.021230.4157D(CNE) does not granger cause SNG1.194950.3149TEL does not granger cause D(CNE)0.868930.5204D(CNE) does not granger cause TEL1.992310.0733TGN does not granger cause D(CNE)1.213910.3051D(CNE) does not granger cause TGN3.895490.0015APX does not granger cause D(D(CNW))0.353650.9063D(D(CNW)) does not granger cause APX0.129340.9924BET-NG does not granger cause D(D(CNW))1.948860.0799D(D(CNW)) does not granger cause BET-NG1.761130.1144EEX-P does not granger cause D(D(CNW))0.435350.8537D(D(CNW)) does not granger cause EEX-P0.663650.6791EEX-B does not granger cause D(D(CNW))0.32820.9208D(D(CNW)) does not granger cause EEX-B0.487020.8167EL does not granger cause D(D(CNW))0.584540.7419D(D(CNW)) does not granger cause EL0.43870.8514SNG does not granger cause D(D(CNW))1.84630.0973D(D(CNW)) does not granger cause SNG1.719580.1238TEL does not granger cause D(D(CNW))1.932760.0824D(D(CNW)) does not granger cause TEL0.969730.4496TGN does not granger cause D(D(CNW))0.767930.5968D(D(CNW)) does not granger cause TGN1.760570.1146Source: author’s own work. Notes: sample: 7/01/2021–12/21/2021. For the definition of variables, please see Table2.The causality test confirms that between the number of new cases of COVID-19 registered in Europe and TEL, respectively TGN, there is a one-way relationship running from CNE to both energy companies listed on the BSE. From an econometric point of view, we came to this decision because thep-value is below the threshold of 10% and 5% respectively, which leads to the rejection of the null hypothesis: the number of new cases of COVID-19 in Europe does not determine Granger-type causation variables analyzed.To our knowledge, the existing studies on the impact of COVID-19 on the energy market during July-December 2021, did not address this type of relationship study, the novelty of the article deepening consisting in this ARDL model that will allow the analysis of long-term relationships between variables selected, as well as the decomposition of the variance and the identification of the causal relations.The innovations of this research include the following aspects. First, this study measured the impact of the COVID-19 pandemic on the energy market. Second, this research integrated advanced econometric models to obtain detailed results on the long-term or short-term relationships between COVID-19 variables and the energy market. Third, this study will help fill the gap in the literature and will be a focal point for future energy market research during the COVID-19 pandemic.


## 4.2.1. Stationarity Analysis

Tables4 and 5 show the outcomes of stationarity analysis for series at the level and in the first difference through the augmented Dickey–Fuller and Kwiatkowski–Phillips–Schmidt–Shin unit root test. A couple of variables are stationary at the first difference (e.g., new cases of COVID-19 worldwide and in Europe) because the probability is above the 1% and 5% relevance level, while some are already stationary (e.g., BET-NG, SNG, EL, TGN, TEL, EEX-B, EEX-P, APX, and LNGI). Alike Bildirici et al. [8], Ali et al. [25], Si et al. [28], Mensi et al. [38], and Yousaf [47], we notice a common integration order of I (0), exception makes the indicators related to the evolution of COVID-19, so in this case, we can reject the null hypothesis and conclude that the series is not stationary.Table 4
The outcomes of the augmented Dickey–Fuller unit root test.
VariablesLevel1st differencet-statisticProb.t-statisticProb.BET-NG−10.968150.00000−8.025040.00000SNG−9.797470.00000−7.022060.00000EL−10.892070.00000−8.212170.00000TGN−11.573740.00000−8.255440.00000TEL−10.629240.00000−7.738480.00000EEX-B−9.284200.00000−10.038590.00000EEX-P−9.262510.00000−7.880230.00000APX−3.852240.00330−9.075060.00000LNGI−10.804120.00000−6.786370.00000D (CNW)−1.054550.73150−6.464690.00000CNE1.492120.99920−3.174390.02410Source: author’s own work. Notes: null hypothesis: each series has a unit root. Intercept included in test equation. Lag length: automatic selection based on Akaike info criterion. Test critical values: 1% level: −3.484198; 5% level: −2.885051; 10% level: −2.579386. For the definition of variables, please see Table2.Table 5
The outcomes of the Kwiatkowski–Phillips–Schmidt–Shin unit root test.
VariablesLevel1st differenceLM-stat.LM-stat.BET-NG0.123170.10802SNG0.051320.27447EL0.073580.22224TGN0.335890.09240TEL0.095700.03197EEX-B0.264770.10671EEX-P0.254600.10521APX0.300140.14376LNGI0.068380.33024D (CNW)0.119720.09233CNE1.172860.13543Source: author’s own work. Notes: null hypothesis: each series is stationary. Intercept included in test equation. Asymptotic critical values: 1% level: 0.739000; 5% level: 0.463000; 10% level: 0.347000. For the definition of variables, please see Table2.Further, Table6 reveals the outcomes of stationarity analysis of the variables with one structural break by means of the Zivot-Andrews unit root test, alike Zhang et al. [7], Le et al. [24]. Accordingly, the occurrence of a structural break in our sample is proved, while the mix integration is strengthened.Table 6
The outcomes of the Zivot–Andrews unit root test.
VariablesLevel1st differencet-statisticProb.Chosen breakpointt-statisticProb.Chosen breakpointBET-NG−11.602500.0145810/29/2021−10.169530.369439/28/2021SNG−10.093280.1099810/28/2021−8.848800.3967611/11/2021EL−11.146200.0322710/28/2021−7.832430.3128411/22/2021TGN−12.346240.0098311/05/2021−9.038700.1122611/24/2021TEL−10.828270.0950610/26/2021−8.099400.0481410/22/2021EEX-B−9.512800.0407710/25/2021−9.625170.1805010/22/2021EEX-P−9.478390.2147810/11/2021−11.075400.1708110/22/2021APX−8.482670.0356010/14/2021−10.477770.0891710/08/2021LNGI−11.552590.0070910/07/2021−8.209490.0262510/04/2021D (CNW)−9.573620.000028/23/2021−15.985630.342039/16/2021CNE−2.455740.0013010/19/2021−6.629530.0239310/06/2021Source: author’s own work. Notes: null hypothesis: each series has a unit root with a structural break in the intercept. Break included in the intercept. 1% critical value: −5.34.5% critical value: −4.93.10% critical value: −4.58. For the definition of variables, please see Table2.


## 4.2.2. Volatility Examination

Because all Jarque–Bera figures are larger than three as Bildirici et al. [8], Zhang and Wang [57] found, and the distributions of the selected series diverge from normally distributed data; GARCH models are required and appropriate for highlighting the progress of volatility. Further, the autocorrelation function (ACF) and partial autocorrelation function (PACF) are used to test for heteroskedasticity. In this regard, Figure 6 exhibits the plots of ACF and PACF. Thus, the existence of serial correlation is inherent in the correlogram of squared returns, except LNGI for which the GARCH model will not be estimated.Figure 6
Plots of autocorrelation function (ACF) and partial autocorrelation function (PACF) for squared returns. Source: author’s own work. Notes: the lag length was set to 20. Variables’ descriptions are provided in Table2.Table7 reports the findings after estimating the GARCH (1, 1) model. In the case of TEL, the sum of ARCH and GARCH parameters is roughly equal to one (e.g., 0.989732), in line with Napon and Asama [4], Bildirici et al. [8]. Hence, at time t, a shock wave will last longer. Specifically, the conditional variance is persistent consistent with Maneejuk et al. [44], Iglesias, and Rivera-Alonso [51]. The findings, however, show the existence of a mean-reverting process because the sum of the ARCH and GARCH effects is less than one. However, in the case of EEX-B, EEX-P, and APX, the sum of parameters is larger than one, thus suggesting that the conditional variance process is explosive.Table 7
The outcomes of the GARCH (1, 1) model.
Dependent variable: BET-NGDependent variable: SNGVariableCoefficientStd. errorz-statisticProb.VariableCoefficientStd. errorz-statisticProb.C0.0005120.0007770.6589560.5099C0.0022260.0009642.307890.021Variance equationVariance equationC4.21E − 054.80E − 050.8768020.3806C9.67E − 050.0001810.5333310.5938RESID (−1)2−0.074930.068867−1.0880460.2766RESID (−1)2−0.048840.072398−0.6746020.4999GARCH (−1)0.6563340.4563521.438220.1504GARCH (−1)0.5420730.9169320.5911810.5544T-DIST. DOF5.0615213.1427831.6105220.1073T-DIST. DOF3.4526231.2507472.7604480.0058ARCH + GARCH0.581404ARCH + GARCH0.493233R-squared−0.001889Mean dependent var8.42E − 05R-squared−0.006586Mean dependent var0.001153AdjustedR-squared−0.001889S.D. dependent var0.009884AdjustedR-squared−0.006586S.D. dependent var0.013269S.E. of regression0.009894Akaike info criterion−6.409628S.E. of regression0.013313Akaike info criterion−5.904976Sum squared resid0.01204Schwarz criterion−6.295907Sum squared resid0.0218Schwarz criterion−5.791255Log likelihood402.3969Hannan–Quinn criterion.−6.363432Log likelihood371.1085Hannan–Quinn criterion.−5.85878Durbin-Watson stat1.988924Durbin-Watson stat1.7568Dependent variable: ELDependent variable: TGNVariableCoefficientStd. errorz-statisticProb.VariableCoefficientStd. errorz-statisticProb.C−8.80E − 060.000686−0.0128180.9898C−0.0015810.000705−2.2426210.0249Variance equationVariance equationC0.0003490.0082760.0421860.9664C6.14E − 065.86E − 061.0488260.2943RESID (−1)2−1.78439542.03255−0.0424530.9661RESID (−1)20.0757880.0746651.0150350.3101GARCH (−1)0.9877170.01889952.262890GARCH (−1)0.8743390.1022928.5474640T-DIST. DOF2.0243890.584223.4651160.0005T-DIST. DOF4.3495662.4443851.7794120.0752ARCH + GARCH−0.796678ARCH + GARCH0.950127R-squared−0.025806Mean dependent var−0.001991R-squared−0.004356Mean dependent var−0.002206AdjustedR-squared−0.025806S.D. dependent var0.012388AdjustedR-squared−0.004356S.D. dependent var0.009498S.E. of regression0.012547Akaike info criterion−6.171115S.E. of regression0.009519Akaike info criterion−6.559459Sum squared resid0.019363Schwarz criterion−6.057394Sum squared resid0.011145Schwarz criterion−6.445738Log likelihood387.6091Hannan–Quinn criterion.−6.124919Log likelihood411.6864Hannan–Quinn criterion.−6.513262Durbin–Watson stat1.93045Durbin–Watson stat2.091561Dependent variable: TELDependent variable: EEX-BVariableCoefficientStd. errorz-statisticProb.VariableCoefficientStd. errorz-statisticProb.C−0.0008560.000926−0.9240740.3554C0.0061580.0107720.5716660.5675Variance equationVariance equationC5.27E − 061.80E − 062.9369290.0033C0.0060220.0147580.4080330.6832RESID (−1)2−0.050780.016305−3.1144110.0018RESID (−1)21.519244.012510.3786260.705GARCH (−1)1.0405120.0139374.696550GARCH (−1)0.671250.1101996.0912320T-DIST. DOF5.7675783.2191841.7916270.0732T-DIST. DOF2.2005450.5998483.6685030.0002ARCH + GARCH0.989732ARCH + GARCH2.19049R-squared−0.000593Mean dependent var−0.001222R-squared−0.000556Mean dependent var0.012439AdjustedR-squared−0.000593S.D. dependent var0.015068AdjustedR-squared−0.000556S.D. dependent var0.267469S.E. of regression0.015073Akaike info criterion−5.758566S.E. of regression0.267543Akaike info criterion−0.16669Sum squared resid0.027944Schwarz criterion−5.644845Sum squared resid8.804247Schwarz criterion−0.052969Log likelihood362.0311Hannan–Quinn criterion.−5.71237Log likelihood15.3348Hannan–Quinn criterion.−0.120494Durbin-Watson stat1.930154Durbin-Watson stat2.342452Dependent variable: EEX-PDependent variable: APXVariableCoefficientStd. errorz-statisticProb.VariableCoefficientStd. errorz-statisticProb.C0.0059190.0116520.5079940.6115C0.0063740.0102060.6245260.5323Variance equationVariance equationC4.012945739.4960.0006990.9994C0.0115820.0459430.2520910.801RESID (−1)2492.9499705091.90.0006990.9994RESID (−1)20.6417692.5009380.2566110.7975GARCH (−1)0.7586570.0865418.7664790GARCH (−1)0.7877820.115036.8484780T-DIST. DOF2.0002990.4277984.6758040T-DIST. DOF2.1486450.6431483.3408230.0008ARCH + GARCH493.708557ARCH + GARCH1.429551R-squared−0.000496Mean dependent var0.013354R-squared−0.001009Mean dependent var0.012429AdjustedR-squared−0.000496S.D. dependent var0.335083AdjustedR-squared−0.001009S.D. dependent var0.191365S.E. of regression0.335166Akaike info criterion0.157424S.E. of regression0.191462Akaike info criterion−0.736736Sum squared resid13.81739Schwarz criterion0.271145Sum squared resid4.508874Schwarz criterion−0.623015Log likelihood−4.760306Hannan–Quinn criterion.0.20362Log likelihood50.67763Hannan–Quinn criterion.−0.69054Durbin–Watson stat2.542814Durbin–Watson stat2.697083Source: author’s own work. Notes: method: ML ARCH–student’st distribution (BFGS/Marquardt steps). Sample: 7/01/2021–12/21/2021. Included observations: 124. For the definition of variables, please see Table 2.We estimated the conditional volatilities (CV) after assessing the GARCH model (1, 1), as shown in Figure7. Among the Romanian companies operating in the energy field, we notice that the highest CV was registered by the major performers in the energy allocation and supply market in Romania, respectively Electrica.Figure 7
Conditional volatility (CV) plot of Romanian stock market index and energy companies. Source: author’s own work. Notes: variables’ descriptions are provided in Table2.


## 4.2.3. ARDL Estimation Results

Since the stationarity test confirmed that the selected variables are integrated in the order of 0 and 1, in line with Le et al. [24], this fact allows us to consider the approach of the ARDL analysis technique. This method permits the cointegration examination of variables that are stationary and non-stationary. For the proper choice of the ARDL model that would allow us to research the relationships that are established between the variables, it is imperative to choose the correct number of offsets. Therefore, alike Nyga-Łukaszewska and Aruga [39], we will analyze the Akaike Information Criteria (AIC) to select the optimal offsets for the variables included in the ARDL model. We will apply the criteria graph, which will indicate the right lags for the ARDL model, and the lowest value is preferred. Figure 8 shows the plots of criteria graph for each ARDL model which considers every energy variable (e.g., BET-NG, SNG, EL, TGN, TEL, EEX-B, EEX-P, APX, and LNGI), as well as the number of new cases at European and global level due to COVID-19. The horizontal axis of each chart represents the ARDL models estimated, and the vertical axis shows the AIC value of the models. The top 20 results are presented in the criteria graph.Figure 8
Akaike information criteria (AIC). Source: author’s own work. Notes: for the definition of variables, please see Table2.
(a)(b)The figures in Table8 signify the results for the ARDL bound test for cointegration. There are provided two critical values for the cointegration test: the lower critical bound assumes all the variables are I (0), meaning that there is no cointegration and the upper bound assumes that all the variables are I (1), meaning that there is cointegration among the variables. If the value of the F-statistic is more than the critical value of bounds, it indicates the long-run relationship between variables. Accordingly, in all cases, the value of the F-statistical test is higher than the limit of 1%, which suggests that there is a long-term relationship between variables. Hence, the null hypothesis is rejected, which means that the variables in all estimated models are cointegrated. As such, we notice cointegration relationships between COVID-19 and BET-NG, SNG, EL, TGN, TEL, EEX-B, EEX-P, APX, and LNGI.Table 8
The results of the ARDL bounds test for the model environment and COVID-19.
Test statistic:F-statisticBET-NG40.32473SNG33.30209EL39.27895TGN52.79222TEL40.24938EEX-B25.06441EEX-P25.75187APX23.92644LNGI38.06681Critical value boundsSignificanceI0 boundI1 bound10%3.174.145%3.794.852.50%4.415.521%5.156.36Source: author’s own work. Notes: null hypothesis: no long-run relationships exist. For the definition of variables, please see Table2.Table9 shows the results of the long-term relationship between variables for COVID-19 new cases in Europe and worldwide and energy market-specific variables. The quantitative outcomes provide support for no impact of the number of new cases of COVID-19 in Europe and globally on the price of energy, natural gas, the BET-NG index, and most of the companies operating in the energy field in Romania. However, there is noticed a positive long-term impact of the number of new cases in the EU and globally on TGN.Table 9
ARDL long-term term coefficients.
VariablesCoefficientStd. errort-statisticProb.APXCNE000.8006750.425D (CNW)00−0.226990.8208BET-NGCNE00−0.997160.3207D (CNW)000.3142190.7539EEX-BCNE000.6295790.5303D (CNW)00−0.438450.6619EEX-PCNE000.6713890.5034D (CNW)00−1.317880.1903ELCNE00−0.825070.411D (CNW)000.8276370.4095LNGICNE000.1109850.9118D (CNW)000.0123880.9901SNGCNE00−0.734280.4642D (CNW)001.5042020.1352TELCNE00−1.351050.1794D (CNW)000.2212260.8253TGNCNE00−2.464090.0152D (CNW)00−2.273020.0249Source: author’s own work. Notes: for the definition of variables, please see Table2.Table10 shows the short-run impact of pandemic on the energy market. For all energy market variables is registered a coefficient of the error correction term (CointEq(−1)) which is negative and significant at the 5% level of significance. Consequently, the negative and significant error correction term, which indicates the speed of conversion, exhibits that on the next day, the dependent variable will reach equilibrium with a speed of between 87% and 229%. Also, the short-term results provided in Table 10 show a positive impact of the new cases of COVID-19 registered in Europe on the evolution of the TEL share price. Therefore, an increase in the number of new cases of COVID-19 in Europe during the period under review leads to an increase in the price of TEL. Such a relationship was also identified in the case of TGN, where the new number of COVID-19 cases both in Europe and globally exerts a positive impact on the share price. Another outstanding result identified from these ARDL models is that, in the short term, the outcomes show a negative impact of new cases of COVID-19 disease in Europe on the variable EEX-B–proxy variable for the price of electricity.Table 10
ARDL cointegrating and short-term coefficients.
VariablesCoefficientStd. errort-statisticProb.Dependent variable: APXD(APX(−1))0.8090480.2277743.5519770.0006D(APX(−2))0.4349990.1630332.6681640.0087D(APX(−3))0.2530360.0916282.7615630.0067D(CNE)000.7959410.4277D(D(CNW))00−0.227210.8207CointEq(−1)−2.2976430.272947−8.417920Dependent variable: BET-NGD(CNE)00−0.993770.3224D(D(CNW))000.3146070.7536CointEq(−1)−1.0037230.09165−10.95170Dependent variable: EEX-BD(EEX-B(−1))0.8332630.1900634.3841530D(EEX-B(−2))0.4648220.1403763.3112570.0013D(EEX-B(−3))0.3138670.0876143.5823950.0005D(CNE)−0.0000010.000001−2.191950.0305D(D(CNW))00−0.439390.6612CointEq(−1)−2.0338930.235973−8.619180Dependent variable: EEX-PD(EEX-P(−1))0.8392160.2083584.0277710.0001D(EEX-P(−2))0.5149260.1524793.3770250.001D(EEX-P(−3))0.2798960.0895423.1258750.0023D(CNE)−0.0000020.000001−1.750910.0827D(D(CNW))00−0.757570.4503CointEq(−1)−2.2064870.257442−8.570810Dependent variable: ELD(CNE)00−0.824380.4114D(D(CNW))000.8330040.4065CointEq(−1)−0.990390.091465−10.82810Dependent variable: LNGID(CNE)000.1109690.9118D(D(CNW))000.0123880.9901CointEq(−1)−0.9831640.091917−10.69620Dependent variable: SNGD(CNE)00−0.737230.4625D(D(CNW))000.8354340.4052CointEq(−1)−0.8734530.089253−9.786280Dependent variable: TELD(CNE)001.3832630.1693D(CNE(−1))000.1019940.9189D(CNE(−2))002.335510.0213D(D(CNW))000.2212640.8253CointEq(−1)−1.0057320.091574−10.98270Dependent variable: TGND(CNE)002.8408560.0053D(CNE(−1))003.2864660.0013D(D(CNW))00−2.363280.0198CointEq(−1)−1.0862950.088083−12.33260Source: author’s own work. Notes: for the definition of variables, please see Table2.Further, alike Jeris and Nath [37], Geyikçi [41], this study applied cumulative sum (CUSUM) and cumulative sum of the squares (CUSUM of Squares) to confirm the stability of the long-run and short-run parameters, respectively the reliability and stability of the examined models. While the cumulative sum of squares test detects abrupt changes from the constancy of the regression coefficient, the cumulative sum test captures systematic variations in the regression coefficients [74]. Under the null hypothesis, the regression coefficients remain constant over time, being equal (or stable) in all sequential subsamples Ploberger and Krämer [75]. Hence, the estimated models are stable and there is no structural break in the observed time series, if the null hypothesis is valid [76]. Contrarily, the alternative assumption is that the regression coefficients fluctuate during the course of the sample. Nevertheless, the test is commonly questioned on the grounds that it is quite straightforward and does not request a priori information regarding the timing of the structural change [78]. Concretely, these tests produce a diagram by recurrently computing the regression coefficients and residuals [76]. For the stability of the estimated outcomes, the graph should maintain the critical values.The results of CUSUM and CUSUM of squares are plotted in Figure9. As far as the CUSUM test is concerned, all the models are stable as the CUSUM blue line is within the 5% significance boundaries illustrated by red dotted lines. Howbeit, the CUSUM test is less powerful because the confidence interval of the test is approximated [77]. In addition, the CUSUM of Squares test shows the stability of APX, BET-NG, EEX-P, and TGN models, but EEX-B, EL, LNGI, SNG, and TEL models are unstable as the CUSUM of Squares blue line crosses the 5% boundary. However, although it is known that the CUSUM of Squares test statistic is distributed as a beta random variable, there should be acknowledged that the confidence interval of the test is also approximated [77].Figure 9
CUSUM test and CUSUM of square test of each ARDL model. Source: author’s own work. Notes: for the definition of variables, please see Table2.
(a)(b)(c)Further, even though several models are unstable as pointed out by Figure9, Table 11 provides supplementary diagnostic tests in order to check for robustness. To avoid the serial correlation problem (the association among included variables and its lagged value), the Breusch-Godfrey LM test for autocorrelation (presence of autocorrelation in the null hypothesis) was conducted, while heteroscedasticity was examined by Breusch-Pagan-Godfrey test. Accordingly, the hypothesis of no serial correlation between variables and its lagged value is rejected at a 5% level of significance (Prob. Chi-Square > 0.05), except for the EEX-P model. In addition, the probability associated with the Chi-Square value is above the 0.05 significance level signifying that the assumption of homoscedasticity fails to be rejected, except EEX-B, EEX-P, and SNG models. Therefore, the errors of the EEX-B, EEX-P, and SNG models are not homoscedastic.Table 11
ARDL diagnosis tests.
APXBreusch–Godfrey serial correlation LM testF-statistic1.956852Prob.F(2, 111)0.1461Obs∗R-squared4.086932Prob. Chi-square(2)0.1296Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic0.871964Prob.F(6, 113)0.5179Obs∗R-squared5.310026Prob. Chi-square(6)0.5047Scaled explained SS13.70353Prob. Chi-square(6)0.0331BET-NGBreusch–Godfrey serial correlation LM testF-statistic0.095708Prob.F(2, 117)0.9088Obs∗R-squared0.200903Prob. Chi-square(2)0.9044Heteroskedasticity test: Breusch-Pagan-GodfreyF-statistic2.362356Prob.F(3, 119)0.0747Obs∗R-squared6.91355Prob. Chi-square(3)0.0747Scaled explained SS10.20941Prob. Chi-square(3)0.0169EEX-BBreusch–Godfrey serial correlation LM testF-statistic1.43912Prob.F(2, 110)0.2416Obs∗R-squared3.059835Prob. Chi-square(2)0.2166Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic3.030858Prob.F(7, 112)0.0059Obs∗R-squared19.11122Prob. Chi-square(7)0.0078Scaled explained SS36.21409Prob. Chi-square(7)0EEX-PBreusch–Godfrey serial correlation LM testF-statistic3.67646Prob.F(2, 109)0.0285Obs∗R-squared7.583398Prob. Chi-square(2)0.0226Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic2.361501Prob.F(8, 111)0.0219Obs∗R-squared17.45327Prob. Chi-square(8)0.0257Scaled explained SS36.76958Prob. Chi-square(8)0ELBreusch–Godfrey serial correlation LM testF-statistic0.078257Prob.F(2, 117)0.9248Obs∗R-squared0.164321Prob. Chi-square(2)0.9211Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic1.515814Prob.F(3, 119)0.214Obs∗R-squared4.527293Prob. Chi-square(3)0.2099Scaled explained SS8.60173Prob. Chi-square(3)0.0351LNGIBreusch–Godfrey serial correlation LM testF-statistic0.03895Prob.F(2, 117)0.9618Obs∗R-squared0.081841Prob. Chi-square(2)0.9599Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic0.019132Prob.F(3, 119)0.9964Obs∗R-squared0.059298Prob. Chi-square(3)0.9962Scaled explained SS1.165012Prob. Chi-square(3)0.7614SNGBreusch–Godfrey serial correlation LM testF-statistic0.210536Prob.F(2, 115)0.8105Obs∗R-squared0.445074Prob. Chi-square(2)0.8005Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic4.599008Prob.F(4, 117)0.0017Obs∗R-squared16.57594Prob. Chi-square(4)0.0023Scaled explained SS38.89259Prob. Chi-square(4)0TELBreusch–Godfrey serial correlation LM testF-statistic0.4525Prob.F(2, 112)0.6372Obs∗R-squared0.969886Prob. Chi-square(2)0.6157Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic0.837925Prob.F(6, 114)0.5431Obs∗R-squared5.110863Prob. Chi-square(6)0.5297Scaled explained SS14.91187Prob. Chi-square(6)0.021TGNBreusch–Godfrey serial correlation LM testF-statistic0.246637Prob.F(2, 114)0.7818Obs∗R-squared0.525616Prob. Chi-square(2)0.7689Heteroskedasticity test: Breusch–Pagan-GodfreyF-statistic1.580195Prob.F(5, 116)0.1711Obs∗R-squared7.779751Prob. Chi-square(5)0.1688Scaled explained SS10.33263Prob. Chi-square(5)0.0663Source: author’s own work. Notes: for the definition of variables, please see Table2.


## 4.2.4. Variance Decomposition Research

To determine the extent to which COVID-19 variables contribute to the explanation of energy market variables, the variance decomposition approach was used, the outcomes being revealed in Table12. Decomposition variation indicates the extent to which a certain variable can explain the evolution of the variation of another variable. In addition, it shows which of the independent variables is stronger in explaining the variability of dependent variables over time.Table 12
The results of the variance decomposition.
Variance decomposition of BET-NGPeriodBET-NGCNED(CNW)498.297860.5159061.186235Variance decomposition of APXPeriodAPXCNED(CNW)495.799620.5795113.620871Variance decomposition of EEX-BPeriodEEX-BCNED(CNW)498.267531.3516390.380828Variance decomposition of EEX-PPeriodEEX-PCNED(CNW)497.445531.8299810.724489Variance decomposition of LNGIPeriodLNGICNED(CNW)497.07622.1712360.75256Variance decomposition of ELPeriodELCNED(CNW)497.454331.6309740.914701Variance decomposition of SNGPeriodSNGCNED(CNW)496.267662.227231.505114Variance decomposition of TELPeriodTELCNED(CNW)498.295610.7208730.983517Variance decomposition of TGNPeriodTGNCNED(CNW)491.53255.3889153.078584Source: author’s own work. Notes: null hypothesis: no long-run relationships exist. Notes: for the definition of variables, please see Table2.From the results obtained after the decomposition of the variance, it can be noticed that the number of new cases of COVID-19 registered in Europe shows a higher contribution in the case of EEX-B, EEX-P, LNGI, EL, SNG, and TGN, whereas for the variables BET-NG, APX, TEL, the largest contribution comes from the number of new COVID-19 cases registered globally.


## 4.2.5. Causality Assessment

To explore causality between selected variables, the Granger causality test is applied. To employ the Granger causality test, the data series must be stationary and have therefore been converted into stationary series. Table13 shows the results after the Granger causality test for the energy market and COVID-19 variables.Table 13
Results of granger causality test.
Null hypothesisF-statisticProb.D(D(CNW)) does not granger cause D(CNE)0.715270.6381D(CNE) does not granger cause D(D(CNW))1.34730.2432APX does not granger cause D(CNE)1.0480.3989D(CNE) does not granger cause APX0.556420.764BET-NG does not granger cause D(CNE)1.013430.4207D(CNE) does not granger cause BET-NG1.564390.1649EEX-P does not granger cause D(CNE)1.561610.1658D(CNE) does not granger cause EEX-P1.454460.2012EEX-B does not granger cause D(CNE)1.856280.0954D(CNE) does not granger cause EEX-B1.261170.2818EL does not granger cause D(CNE)1.020950.4159D(CNE) does not granger cause EL0.858160.5284SNG does not granger cause D(CNE)1.021230.4157D(CNE) does not granger cause SNG1.194950.3149TEL does not granger cause D(CNE)0.868930.5204D(CNE) does not granger cause TEL1.992310.0733TGN does not granger cause D(CNE)1.213910.3051D(CNE) does not granger cause TGN3.895490.0015APX does not granger cause D(D(CNW))0.353650.9063D(D(CNW)) does not granger cause APX0.129340.9924BET-NG does not granger cause D(D(CNW))1.948860.0799D(D(CNW)) does not granger cause BET-NG1.761130.1144EEX-P does not granger cause D(D(CNW))0.435350.8537D(D(CNW)) does not granger cause EEX-P0.663650.6791EEX-B does not granger cause D(D(CNW))0.32820.9208D(D(CNW)) does not granger cause EEX-B0.487020.8167EL does not granger cause D(D(CNW))0.584540.7419D(D(CNW)) does not granger cause EL0.43870.8514SNG does not granger cause D(D(CNW))1.84630.0973D(D(CNW)) does not granger cause SNG1.719580.1238TEL does not granger cause D(D(CNW))1.932760.0824D(D(CNW)) does not granger cause TEL0.969730.4496TGN does not granger cause D(D(CNW))0.767930.5968D(D(CNW)) does not granger cause TGN1.760570.1146Source: author’s own work. Notes: sample: 7/01/2021–12/21/2021. For the definition of variables, please see Table2.The causality test confirms that between the number of new cases of COVID-19 registered in Europe and TEL, respectively TGN, there is a one-way relationship running from CNE to both energy companies listed on the BSE. From an econometric point of view, we came to this decision because thep-value is below the threshold of 10% and 5% respectively, which leads to the rejection of the null hypothesis: the number of new cases of COVID-19 in Europe does not determine Granger-type causation variables analyzed.To our knowledge, the existing studies on the impact of COVID-19 on the energy market during July-December 2021, did not address this type of relationship study, the novelty of the article deepening consisting in this ARDL model that will allow the analysis of long-term relationships between variables selected, as well as the decomposition of the variance and the identification of the causal relations.The innovations of this research include the following aspects. First, this study measured the impact of the COVID-19 pandemic on the energy market. Second, this research integrated advanced econometric models to obtain detailed results on the long-term or short-term relationships between COVID-19 variables and the energy market. Third, this study will help fill the gap in the literature and will be a focal point for future energy market research during the COVID-19 pandemic.


## 5. Concluding Remarks and Policy Implications

The extreme insecurity of the pandemic and the related economic failures have made markets very volatile and unpredictable. Hence, the risks of the global financial market boosted considerably in response to the disease. In this article, we examined how the number of new COVID-19 cases in Europe and globally is affecting the Romanian energy market. To achieve our goal, we used daily data for the period July 1, 2021–December 21, 2021, which includes the fourth wave of the COVID-19 pandemic. We have selected a wide range of variables that characterize the energy market: energy price indices, natural gas, and the BET-NG index of the Bucharest Stock Exchange. Also, we covered several companies listed on BSE that act in the energy field which are also among the most traded companies in Romania.Due to the mix of stationary and nonstationary variables, the ARDL model was adopted in this research. Through this model, it was feasible to study the relationships that are established in the long and short term. There should be noted that the number of short-term relationships was higher, which is also due to the short episodes of high volatility recorded in the investigated timeframe. However, in the long run, there is a positive impact on the number of new cases in Europe and in the world on the return of TGN. In contrast, no long-term relationship was identified between COVID-19 variables and the price of electricity and natural gas. In the short term, a negative impact of new cases of COVID-19 infection in Europe on the price of electricity has been identified. As well, in the short run, the positive impact of the new COVID-19 cases is noticed in the companies traded on BSE. Hence, an increase in the number of new COVID-19 cases in Europe during the analyzed period leads to a rise in the share price of TEL. In the same vein, an increase in the number of new cases of COVID-19 both in Europe and globally have a positive impact on the price of TGN shares.To explore causality between variables, the Granger causality test was applied. One-way causal relationships have been identified from the number of new COVID-19 cases in Europe to the returns of TEL and TGN. Nevertheless, no Granger causal relationships have been found between the COVID-19 variables and the price of electricity or natural gas.This research showed that the variables do not have a direct impact on energy prices, but certainly, the effects of the COVID-19 pandemic are those that contributed to the increase in energy prices, obviously indirectly by reflecting the negative effects of measures to combat the spread of COVID-19 virus. Our empirical findings provide insight into how the energy market is affected by the COVID-19 pandemic.As it turned out, prices have risen, and this is due to the supply deficit. Therefore, we consider a series of measures that should be implemented by the administrative sector in order to avoid a social crisis. In the immediate vicinity, aid may be provided to vulnerable household consumers for price increases, either in cash or by partial payment of the budget bill; temporary delays in paying bills and avoiding disconnections during this period; reduction of taxes and duties (for the most vulnerable consumers); state aid to industry; increasing market liquidity and transparency; and investigations into speculation or market abuse.We also presume that increasing investment in renewables, as well as energy efficiency in buildings and industry, is a solution that will only bring positive results in the long run. As a result, transitioning to clean energy is the best way to protect against future price shocks and should be accelerated. Therefore, the findings of this study can help investors optimize their portfolios while also providing guidance to decision-makers and regulators. In addition, in order to help the oil market, governments may relax quarantine restrictions and reopen their businesses [12].Because of the ongoing changes in the pandemic, as well as the various geopolitical and economic events that may occur, more research on this topic is likely in the future. These findings should be considered by investors and policymakers because they argue that the relationship between energy goods and COVID-19 is dynamic rather than linear. For the reason that the current study is limited to the fourth wave of the disease, upcoming investigations should also cover the preceding COVID-19 waves. The current study is limited to merely the Romanian energy market and the fourth wave of the COVID-19 pandemic. However, further research avenues should cover more international energy markets, as well as the prior pandemic period. As well, upcoming studies might assess the co-movements of the energy market with other variables such as gold or agricultural commodities.



---

*Source: 3380183-2022-09-26.xml*