Modelling volatility in job loss during the COVID-19 pandemic: The Malaysian case

Abstract

This study employs a suitable volatility model that examines the impact of COVID-19 new cases and deaths on the volatility of daily job loss in Malaysia. Autoregressive Distributed Lag (ARDL) and Generalized Autoregressive Conditional Heteroscedasticity (GARCH) were employed as the modelling strategy to estimate daily data from January to December 2020. In addition, the asymmetric GARCH-M (EGARCH-M, TGARCH-M, and PGARCH-M) were further applied. The findings from different versions of the ARDL(p,q1,q2)-(E,T,P)GARCH(1,1)-M model show that the ARDL-EGARCH-M model can capture the volatility and clustering of variability in job loss. The findings revealed asymmetry effects, suggesting that negative shocks (bad news) in a pandemic period increased volatility in job loss compared to positive shocks (good news). Policy implications relating to lockdown measures and news signals were provided.

Keywords:

• Covid-19
• job loss
• ARDL
• GARCH-M
• Malaysia

JEL Classification:

• E20
• I10
• J64

1. Introduction

The world has seen the devastating impact of the COVID-19 pandemic by crippling most economic activities. The unprecedented pandemic has not only caused a global health disaster but it has also caused an economic and labour market crisis. To counter this pandemic from spreading, countries around the world have implemented measures that slow the spread of the coronavirus. At the same time, measures, such as the non-pharmaceutical containment has constraint economic growth (World Bank, 2020). According to Cotofan et al. (2021), the global economic growth is expected to contract by more than 4% in 2020, and the most affected sectors of the economy include tourism and hospitality related sectors such as hotels, restaurants, wholesales and retails, crafts and shopping malls, movie theatres, cafes, airlines, and other land and sea modes of transportation.

The containment measures enforced by the governments around the world have had a significant impact on the worldwide labour market (Cotofan et al., 2021). Business closures, loss of employment, higher unemployment rates, lower labour participation rates, and reduced hours worked (OECD, 2021) were some of the negative effects of the COVID-19 pandemic (Almeida & Santos, 2020; Deady et al., 2020; Fontaine, 2020; Kong & Prinz, 2020). According to Van der Wielen and Barrios (2020), the COVID-19 pandemic intensified people’s economic anxiety, resulting in a significant slowdown in European Union nations’ labor markets and spending. The COVID-19 pandemic has utterly shattered any sense of security at work, and the labor market sector is presently witnessing unprecedented job losses (International Labor Organization, 2020).

Similar events in different countries had varying consequences on the job market. Female, young, low-paid workers, and members of various ethnic minority groups were among those who lost their jobs in the United Kingdom as a result of the lockdown imposed by the UK government to prevent the COVID-19 virus from spreading (Blundell et al., 2020). According to Couch et al. (2020), the unemployment rate for African-Americans in the United States (US) grew to 16.6%, while it was 18.2% for Latinx. They argued that because of adverse occupational distribution and weaker talents, Latinx have had much higher unemployment rates than Whites. Furthermore, Béland et al. (2020) revealed that the COVID-19 had a more severe negative impact on the US labor market. According to Falk et al. (2020), the COVID-19 pandemic led to high unemployment rates in the United States for young workers, women, workers with low educational attainment, part-time workers, and racial and ethnic minorities.

Malaysia has not been spared the COVID-19 shock. The threat projected by COVID-19 became a reality on 25 January 2020, when the first COVID-19 positive cases were revealed. It has been claimed that COVID-19 and the Malaysian government’s lockdown tactics have resulted in population job losses (Shah et al., 2020). The Malaysian government began implementing its mobility control order (MCO), or “lockdown” measures, on 18 March 2020. Shutting non-essential businesses, closing schools and workplaces, issuing a stay-at-home order, prohibiting large-scale gatherings and public events, and restricting domestic and foreign travel are among the measures. According to the Department of Statistics Malaysia (DOSM), the unemployment rate climbed from 3.3% in February to 3.9% in March 2020, up from 3.3% in February and 3.2% in January 2020 (Department of Statistics Malaysia, 2020a). In May, the jobless rate peaked at 5.3%. It subsequently began to decline, hitting 4.7% in August 2020 before slightly climbing to 4.8% in December 2020.

Nonetheless, COVID-19 has had a disproportionately detrimental impact on the Malaysian labour market, as it has in many other countries. The employment rate has fallen from 69.1% in the fourth quarter of 2019 to 68.1% in the second quarter of 2020. There were more female jobless (5.5%) than male jobless (4.7%) in the second quarter of 2020, as well as more youthful employees (15–24 years old) than senior workers (25–34 years old) (Department of Statistics Malaysia, 2020c). The unemployment rate for men and women is equal in the fourth quarter of 2020, at 4.8% each, but it climbed to 12.7% for young employees (15–24 years old) and 4.8% for older workers (25–34 years old) (Department of Statistics Malaysia, 2020b).

Figures 1 and 2 provide time plots for the number of people unemployed, new COVID-19 cases, and new deaths. Figure 1 depicts the number of people who lost their job on a daily basis. The number of jobs lost during the first COVID-19 wave peaked on 2 July 2020, when 1,540 were unemployed. We can clearly observe that the number of new COVID-19 cases, deaths, and job losses will alter dynamically during the year 2020. Following the Sabah state election on 26 September 2020, we anticipate an increase in the number of new COVID-19 cases that will persist through the end of 2020. The COVID-19 new deaths volatility, on the other hand, can be shown to be focused in two-time frames: the first, from mid-March to mid-May 2020, and the second, from mid-September to the end of 2020. There were 19 further deaths recorded after September. This is directly tied to the increase in new COVID-19 cases following the 2020 Sabah state election.

Figure 1. Daily number of job loss and number of COVID-19 new cases.

Figure 2. Daily number of job loss and number of COVID-19 new deaths.

The paper is inspired with the following reasons, notably: (i) the unprecedented and alarming increase in the number of COVID-19 new cases and death cases; (ii) the negative impact of government lockdown measures on job losses especially among the vulnerable income segments; (iii) the disproportional effect of COVID-19 pandemic and lockdown measures on the Malaysian labor market; (iv) the necessity to identify the asymmetric effect of COVID-19 pandemic on the volatility of job loss; (v) and the gap in extant literature of not modelling volatility in job loss during the COVID-19 era.

This paper potentially contributes to the literature as it is different from recent related studies such as the distributional effect of lockdown in Malaysia (Saari et al., 2022), impact of COVID-19 on SMEs operation (Shafi et al., 2020), lockdown, distribution of income and food security (Arndt et al., 2020), income support on job search, COVID-19 and workplace mobility (Asfaw, 2021), COVID-19 effect on agriculture and food security (Adhikari et al., 2021). The purpose of this research is to find an appropriate model that accounts for the influence of COVID-19 new cases and deaths on the daily volatility of job loss in Malaysia from January to December 2020. To determine how Malaysian labour market respond to the outbreak, we used daily data on loss of employment as a stand-in for daily unemployment rates. Since daily macroeconomic time series are characterised by volatility and clustering, we apply the Generalized Autoregressive Heteroscedasticity (GARCH) to cope with volatility in predicting job loss in Malaysia from January to December 2020.

This paper is divided into five sections. Section 2 follows the Introduction and provides basic statistics for the data used in this work. Section 3 presents estimation methods, including symmetric and asymmetric GARCH models. Section 4 evaluates the models with the purpose of determining the model that best depicts the unpredictability of job loss. Section 5 contains the conclusion.

2. Theoretical background

Originally, the COVID-19 pandemic resulted in a global health crisis before it transformed the economy into a global economic crisis. From an economic perspective, health can affect economic growth by providing healthy people to performed productive activities. Workers with better health will be able to do work better, low absentee, acquire higher education and skills, reduce cost of illnesses and increase productivity (Van der Hout, 2015). Thus, health can be considered as one of the important components of human capital stock that generates output. In fact, Barro (2013) shows that a rise in life expectancy from 50 to 70 years would raise growth rate by 1.4% per year.

Using a closed-economy model with overlapping-generations, Augier and Yaly (2013) present that for a country with less severe disease level, the disease will progressively disappear. But, if the disease is above the threshold level, the disease will spread rapidly and the economy will be vulnerable to facing economic crisis. The event study conducted by Goda and Soltas (2022) show that workers who miss an entire week due COVID-19 illness are approximately 7% points less likely to be in the labour force one year after absence compared to workers who do not miss work for health reasons. They further reiterated that COVID-19 illnesses have reduced labour force participation rates by approximately half a million workers through June 2022.

On the other hand, Xiang et al. (2021) combined economic theory and infectious disease model and analyses the long-term impact of the COVID-19 pandemic on labour supply, household consumption, and economic growth. They found that the outbreak of COVID-19 directly affects labour supply and output, resulting in a significant negative impact on the economic growth. They further expound that only when the government increases investment in the prevention measures and medical treatment at the same time, the negative consequences of pandemic are able to be minimized. Goes and Gallo (2021) provide a theoretical framework based on a two-equations dynamical system that capture the interactions between the COVID-19 and unemployment rate. Their findings suggest that to suppress to effects of the Covid-19, the containment measures in the early phase of the disease spreading would have to be exactly as strong as the force of infection.

The COVID-19 pandemic had significant effects on both the supply and demand of the labour force. Supply and demand shocks are likely intertwined with what starts as a supply shock—lockdown, layoffs and firm exit—leading to demand shock. Guerrieri et al. (2022) showed that, in a two-sector new Keynesian model with low substitutability in consumption, asymmetric labour supply shocks can lead to reductions in demand that are higher than the initial shock. The model further demonstrates that supply shocks can lead to decreases in demand. Del Rio-Chanona et al. (2020) exert that one-fifth of the US economy’s output, one-quarter of consumption, and 16% of the US total wage income will be reduced as a result of the COVID-19 pandemic. They emphasis that supply shocks account for the majority of the reduction as large part of the manufacturing and services sectors being classified as non-essential, while its labour force is unable to work from home.

Cakmakli et al. (2020) point out that countries that have greater backward linkages would likely be affected by the supply shock, while countries that have greater forward linkages will be affected by demand shock. For example, countries that export raw resources and materials for manufacturing are likely to experience large demand shock while Novy and Jary (2021) point out that COVID-19 affect the demand side of the economy. All components of aggregate demand—consumption, investment, government expenditures, and foreign trade are adversely affected. The social distancing and lockdown measures have resulted in the closure of business, reduction in household consumption, income, wealth, job losses, and reduced wages. On the other hand, theoretical model developed by Eichenbaum et al. (2020) demonstrate that reductions in consumption and work will slow down the spread of the Covid-19, but at the same time, the measures will also increase the severity of economic recession.

The lockdown measures have cause job loss and income in both the formal and informal sectors. As reported in The Financial Express (cited in Hossain, 2021), job lost as a result of the COVID-19 pandemic was about 7 million peoples lost their job in the micro, small and medium enterprises (MSME) sector, 5 million jobs in the transport sector, 3 million jobs in the construction sector and 3 million jobs in the manufacturing sector. Bauer and Weber (2021) found that 60% of the increased inflows from employment into unemployment in April 2020 were due to the shutdown measures. The shutdown measures increased unemployment in the short-run by 117,000 people.

However, Su et al. (2021) showed that in five of the selected European countries, COVID-19 new cases cause unemployment in Germany, Italy and the UK, while COVID-19 deaths cause unemployment in Italy and UK. In the African continent, the work by Ranchhod and Daniels (2021) demonstrated that COVID-19 pandemic increased the unemployment rate rapidly. Similarly, in Asia, by employing an advanced hybrid modeling approach, Shi et al. (2022) found that the unemployment rate in seven developing countries of Asia: Iran, Sri Lanka; Bangladesh; Pakistan; Indonesia; China; and India.

3. Methodology

3.1. Modelling unemployment volatility

The present study’s main purpose is to evaluate how the regressors COVID-19 new cases and new deaths interact with job loss. In this study, we intend to use the GARCH model to investigate the association between daily job losses and daily news on new cases and new deaths connected with Covid-19. A GARCH model must include the mean equation and the variance equation, both of which are short-run equations.

Before estimating the mean equation for job loss, we must define the long-term relationship between the variables. Say we compute the subsequent long-run level regression.

(1) $\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}={\psi }_0+{\psi }_1\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}+{\psi }_2\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}}+{\epsilon }_{\mathrm{t}}$(1)

where $\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}$ is the logarithm of job loss, $\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}$ and $\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}}$ are the logarithm of number of new cases, and number of new deaths; and the error term ${\epsilon }_{\mathrm{t}}$ is assumed to be white noise.

Before estimating Equation (1)(1) $\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}={\psi }_0+{\psi }_1\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}+{\psi }_2\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}}+{\epsilon }_{\mathrm{t}}$(1) that contains integrated series, testing for the order of integration of a series is an important task in economic time series. When non-stationary variables are regressed, they provide erroneous regression results. Given this, the standard augmented Dickey–Fuller unit root test devised by Dickey and Fuller (1981) is employed in this study to determine the order in which the series under consideration were integrated. The null hypothesis of a unit root will be evaluated first on the level of the series (including the deterministic term—intercept or/and trend), and then on their first-differences. If the null hypothesis cannot be rejected at the level but the unit root can be rejected in first-difference, we may infer that the series ${\mathrm{y}}_{\mathrm{t}}$ is non-stationary in levels but becomes stationary after first-differences. In other words, ${\mathrm{y}}_{\mathrm{t}}\text{~}\mathrm{I}\left(1\right)$ and $\mathrm{\Delta }{\mathrm{y}}_{\mathrm{t}}\text{~}\mathrm{I}\left(0\right)$.

We can proceed with the cointegration test after determining, for example, that the number of new cases, new deaths, and job losses are all I(1) in level. The main purpose of conducting cointegration for the three variables is to determine the applicability of the long-run model as stated in Equation (1)(1) $\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}={\psi }_0+{\psi }_1\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}+{\psi }_2\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}}+{\epsilon }_{\mathrm{t}}$(1) . The cointegrability of the variables ensures that the estimated regression is not spurious. For a cointegrated model, all statistical properties are valid, and conclusions may be formed using standard statistical indicators.

In this paper, we employ the Autoregressive Distributed Lag (ARDL) approach proposed by Pesaran et al. (2001). When applied to a mixture of I(0) and I(1) variables, small samples, endogeneity, and models with suitable lag structures, the ARDL approach is effective and trustworthy. Furthermore, Pesaran et al. (2001) show that utilising the ARDL technique, simultaneous estimate of the long-run and short-run (our mean equation) models is conceivable. This framework also supports the Bound-F test for cointegration. According to Pesaran et al. (2001), the ARDL(1,1) model in levels may be utilised to build a long-run model as shown in Equation (1)(1) $\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}={\psi }_0+{\psi }_1\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}+{\psi }_2\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}}+{\epsilon }_{\mathrm{t}}$(1) . With the lag operator L, we have $\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}-1}$=$\mathrm{L}\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}$ so that:

(2) $\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}={\psi }_0+{\psi }_1\left(\mathrm{L}\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}\right)+{\psi }_2\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}+{\psi }_3\left(\mathrm{L}\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}\right)+{\psi }_4\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}}+{\psi }_5\left(\mathrm{L}\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}}\right)+{\eta }_{\mathrm{t}}$(2)

Thus, Equation (2)(2) $\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}={\psi }_0+{\psi }_1\left(\mathrm{L}\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}\right)+{\psi }_2\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}+{\psi }_3\left(\mathrm{L}\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}\right)+{\psi }_4\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}}+{\psi }_5\left(\mathrm{L}\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}}\right)+{\eta }_{\mathrm{t}}$(2) can be rewritten as:

$\left(1-\mathrm{L}{\psi }_1\right)\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}={\psi }_0+({\psi }_2+\mathrm{L}{\psi }_3)\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}+({\psi }_4+\mathrm{L}{\psi }_5)\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}}+{\eta }_{\mathrm{t}}$

(3) $\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}={\psi }_0\left(1-\mathrm{L}{\psi }_1\right)+\left[({\psi }_2+\mathrm{L}{\psi }_3\right)/\left(1-\mathrm{L}{\psi }_1\right)\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}+\left[({\psi }_4+\mathrm{L}{\psi }_5\right)/\left(1-\mathrm{L}{\psi }_1\right)\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}}+{\eta }_{\mathrm{t}}$(3)

Equation (2)(2) $\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}={\psi }_0+{\psi }_1\left(\mathrm{L}\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}\right)+{\psi }_2\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}+{\psi }_3\left(\mathrm{L}\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}\right)+{\psi }_4\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}}+{\psi }_5\left(\mathrm{L}\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}}\right)+{\eta }_{\mathrm{t}}$(2) must pass the non-serial correlation test with optimum lag length.

Nonetheless, Pesaran et al. (2001) proposed the Bounds F-test on the following conditional ARDL-error-correction model (ARDL-ECM) to test for cointegration.

$\mathrm{\Delta }\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}={\rho }_0+{\rho }_1\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}-1}+{\rho }_2\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-1}+{\rho }_3\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}-1}\hspace{0.17em}+\sum _{\mathrm{i}=1}^{\mathrm{p}}{\vartheta }_{1\mathrm{i}}\mathrm{\Delta }\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}-\mathrm{i}}$

(4) $+\sum _{\mathrm{i}=0}^{\mathrm{q}}{\vartheta }_{2\mathrm{i}}\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-\mathrm{i}}+\sum _{\mathrm{i}=0}^{\mathrm{q}}{\vartheta }_{3\mathrm{i}}\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}-\mathrm{i}}+{\in }_{\mathrm{t}}$(4)

The bound-F tests were conducted to see if ${\rho }_i=0$ (null hypothesis) or at least one ${\rho }_i$ is not equal to zero (alternative hypothesis). The long-run cointegrating relationship is recognised when the estimated F-statistic is compared to the bound critical values supplied by Narayan (2005) for small sample size. When the computed F-statistic exceeds the upper boundaries of the critical value for cointegration, the null hypothesis of no cointegration is rejected. The variables are not cointegrated; however, if the computed F-statistic is less than the critical value lower limits and the null hypothesis of no cointegration is not rejected.

If the calculated F-statistic falls between the upper and lower bounds of the critical value, the judgement is inconclusive. The non-cointegration null hypothesis was rejected, suggesting that cointegration exists and that Equation (1)(1) $\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}={\psi }_0+{\psi }_1\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}+{\psi }_2\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}}+{\epsilon }_{\mathrm{t}}$(1) is a trustworthy, non-spurious long-run model.

After estimating the long-run cointegrating regression, the short-run model, or error-correction model, can be stated as follows:

$\mathrm{\Delta }\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}={\delta }_0+\pi EC{\mathrm{M}}_{\mathrm{t}-1}+\sum _{\mathrm{i}=1}^{\mathrm{p}}{\delta }_{1\mathrm{i}}\mathrm{\Delta }\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}-\mathrm{i}}+\sum _{\mathrm{i}=0}^{\mathrm{q}}{\delta }_{2\mathrm{i}}\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-\mathrm{i}}$

(5) $+\sum _{\mathrm{i}=0}^{\mathrm{q}}{\delta }_{3\mathrm{i}}\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}-\mathrm{i}}+{\mu }_{\mathrm{t}}$(5)

where $\mathrm{E}\mathrm{C}{\mathrm{M}}_{\mathrm{t}-1}={\epsilon }_{\mathrm{t}-1}=\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}-1}-\left[{\psi }_0+{\psi }_1\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-1}+{\psi }_2\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}-1}\right]$. The significance of the estimated coefficient and negative values would also imply cointegration (Engle & Granger, 1987). The estimated coefficient’s range is between 0 and −2, according to Loayza and Rancière (2006), Blanco (2013), and Fromentina and Leon (2019). The error-correction short-run model is unique in that it includes long-run information on both $\mathrm{j}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{s}{\mathrm{s}}_{\mathrm{t}}$ and $\mathrm{c}\mathrm{o}\mathrm{v}\mathrm{i}{\mathrm{d}}_{\mathrm{j}\mathrm{t}}$ in the short-run model, which is I(0) as represented by the $\mathrm{E}\mathrm{C}{\mathrm{M}}_{\mathrm{t}-1}$term.

After obtaining our mean equation (also known as the error-correction model), we may specify the variance equation. In this paper, we employ the GARCH and exponential GARCH models and their variants, which can account for non-constant variance over time, to account for the volatility of the job loss series, as shown in Figures 3(3), 4 (4), and 5(5), and empirical evidence shown in Table 1. The two types of GARCH models are the symmetric GARCH model and the asymmetric GARCH model.

Figure 3. Daily log and log changes in number of job loss, number of COVID-19 new cases and new deaths.

Figure 4. Q-Q plots of log job loss, COVID-19 new cases and new deaths, and log changes in loss of employment, COVID-19 new cases and new deaths.

Figure 5. Q-Q plot of job loss regression residuals.

Table 1. Descriptive statistics

Series	No. obs	Mean	Max	Min	Std. Dev.	Skewness	Kurtosis	Jarque-Bera	Anderson-Darling
Job loss	363	6.11	8.03	2.64	0.80	−0.84	4.46	75.28***	3.59***
New cases	363	4.30	8.53	0.00	2.55	−0.16	2.02	16.11***	4.46***
New deaths	363	0.69	3.18	0.00	0.92	0.92	2.43	56.30***	41.63***
$\mathrm{\Delta }\mathrm{J}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}$	363	0.00	2.74	−3.09	0.72	0.15	5.14	70.35***	7.02***
$\mathrm{\Delta }\mathrm{N}$ew cases	363	0.02	4.56	−3.69	0.84	0.51	8.98	556.0***	10.53***
$\mathrm{\Delta }\mathrm{N}$ew deaths	363	0.01	2.31	−3.09	0.67	−0.09	6.22	157.2***	29.27***

Notes: Asterisks *** denotes statistically significant at 1% level. Series of job loss, new cases and new deaths are in natural logarithm.

3.2. The symmetric GARCH model

The majority of volatility modelling research has focused on financial time series. Autoregressive conditional heteroscedasticity (ARCH) and its extension are the key methodologies utilised in modelling and predicting stock market volatility (GARCH). The GARCH model, introduced by Bollerslev (1986), may account to explain volatility clustering. Under the GARCH model, the conditional variance can be affected by its own prior delays. The mean and variance equations are two distinct specifications required by each GARCH family model. In most situations, a GARCH(1,1) was sufficient to find volatility clustering in the data (Engle, 2004). The specification of the estimated GARCH(1,1) with conditional mean equation shown in Table 3 are as follows:

$\mathrm{\Delta }\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}={\delta }_0+\pi EC{\mathrm{M}}_{\mathrm{t}-1}+{\delta }_1\mathrm{\Delta }\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}-1}+{\delta }_2\mathrm{\Delta }\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}-2}+{\delta }_3\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}$

(6) $+{\delta }_4\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-1}+{\delta }_5\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-2}+{\in }_{\mathrm{t}},{\in }_{\mathrm{t}}{ }_{{ }_{\text{~}}}\left(0,{\sigma }_{\mathrm{t}}^2\right)$(6)

(7) ${\sigma }_{\mathrm{t}}^2=\omega +\alpha {\text{isin}}_{\mathrm{t}-1}^2+{\beta \sigma }_{\mathrm{t}-1}^2$(7)

However, as shown in Table 3, the error correction model or our mean equation is not free of the heteroscedasticity problem. To overcome this issue, Engle et al. (1987) proposed that the conditional mean be a function of the conditional variance in the following way,

$\mathrm{\Delta }\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}}={\delta }_0+\pi EC{\mathrm{M}}_{\mathrm{t}-1}+{\delta }_1\mathrm{\Delta }\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}-1}+{\delta }_2\mathrm{\Delta }\mathrm{l}\mathrm{o}{\mathrm{e}}_{\mathrm{t}-2}+{\delta }_3\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}$

(8) $+{\delta }_4\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-1}+\hspace{0.17em}{\delta }_5\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-2}+{\lambda \sigma }_{\mathrm{t}}^2+{\text{isin}}_{\mathrm{t}},\hspace{0.17em}{\text{isin}}_{\mathrm{t}}\text{~}\left(0,{\sigma }_{\mathrm{t}}^2\right)$(8)

The mean Equation (8)(8) $+{\delta }_4\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-1}+\hspace{0.17em}{\delta }_5\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-2}+{\lambda \sigma }_{\mathrm{t}}^2+{\text{isin}}_{\mathrm{t}},\hspace{0.17em}{\text{isin}}_{\mathrm{t}}\text{~}\left(0,{\sigma }_{\mathrm{t}}^2\right)$(8) is known as the GARCH-in-Mean (GARCH-M). It illustrates that changes in job loss are influenced by a constant ${\delta }_0$; the linear combination of job loss, COVID-19 new cases, and new deaths (i.e., the long-run model residuals lagged one period), changes in job loss with lagged one and two periods, current changes in COVID-19 new cases, changes in new cases with lagged one and two periods, and the conditional variance (volatility or shock). The variance Equation (7)(7) ${\sigma }_{\mathrm{t}}^2=\omega +\alpha {\text{isin}}_{\mathrm{t}-1}^2+{\beta \sigma }_{\mathrm{t}-1}^2$(7) , on the other hand, implies that the conditional variance of ${\sigma }_{\mathrm{t}}^2$ is reliant on both its conditional variance lagged one period (${\sigma }_{\mathrm{t}-1}^2$) and its squared error delayed one period $\left({\text{isin}}_{\mathrm{t}-1}^2\right)$.

Long-term average volatility is a constant $\omega$; current volatility news $\left(\alpha \right)$ and historical volatility information $\left(\beta \right)$, on the other hand, show how volatility is impacted. To ensure that volatility is always positive, the parameters $\omega$, $\alpha$, $\beta$ and $\lambda$ are assumed to be non-negative. Furthermore, the stationary condition for GARCH(1,1) is $\alpha +\beta <1$, and as when $\alpha +\beta$ approaches 1, the shock to volatility decays at a slower pace. When the sum of the ARCH and GARCH terms is close to one, the volatility is persistent, which means it may take longer for the volatility to revert to a quieter phase. Furthermore, the unconditional variance of ${\text{isin}}_{\mathrm{t}}$ is undefinable for $\alpha +\beta >1$, a situation known as “non-stationarity in variance,” whereas $\alpha +\beta =1$ is known as a “unit root in variance.” A positive and significant, on the other hand, indicates that more volatility, as represented by ${\sigma }_{\mathrm{t}}^2$, leads to an increase in changes in job loss or vice versa.

3.3. The asymmetric GARCH models

3.3.1. The EGARCH model

The conditional variance cannot respond asymmetrically to changes in the volatile series, which is a limitation of the GARCH model. According to leverage effects, the conditional variance ${\sigma }_{\mathrm{t}}^2$ might react asymmetrically to both positive and negative values of $\mathrm{j}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{s}{\mathrm{s}}_{\mathrm{t}}$. Nelson (1991) proposes the EGARCH model as a solution to avoid the symmetrical GARCH by representing asymmetric reactions of the time-varying variance to shocks while simultaneously ensuring that the variance is always positive. An EGARCH(1,1) model is defined as follows.

(9) $log\left({\sigma }_{\mathrm{t}}^2\right)=\omega +\alpha \left|\frac{{\text{isin}}_{\mathrm{t}-1}}{\sqrt{{\sigma }_{\mathrm{t}-1}^2}}\right|+\gamma \frac{{\text{isin}}_{\mathrm{t}-1}}{\sqrt{{\sigma }_{\mathrm{t}-1}^2}}+\beta log\left({\sigma }_{\mathrm{t}-1}^2\right)$(9)

where the conditional variance’s logarithm is on the left side of Equation (9)(9) $log\left({\sigma }_{\mathrm{t}}^2\right)=\omega +\alpha \left|\frac{{\text{isin}}_{\mathrm{t}-1}}{\sqrt{{\sigma }_{\mathrm{t}-1}^2}}\right|+\gamma \frac{{\text{isin}}_{\mathrm{t}-1}}{\sqrt{{\sigma }_{\mathrm{t}-1}^2}}+\beta log\left({\sigma }_{\mathrm{t}-1}^2\right)$(9) . As a result, EGARCH imposes no non-negative restrictions on the model parameters $\omega$, $\alpha$, $\gamma$ and $\beta$. This suggests that the leverage effect is exponential rather than quadratic, and hence conditional variance forecasts are guaranteed to be non-negative. However, in order to maintain the covariance stationary, $\beta$ must be positive and less than 1. While $\beta$ measures the persistence of conditional volatility, the parameters $\alpha$ and $\gamma$ characterise the model’s magnitude (or size) effect and symmetric effect, known as the GARCH effects. This implies that volatility continues for a long period following a “market crisis” when $\beta$ is significantly large (Alexander, 2009).

The leverage impact or the asymmetric repercussions of volatility shocks are measured by this parameter $\gamma$. When $\gamma$ = 0, the model is symmetric; when $\gamma \ne 0$, the leverage effect is present. A value of zero indicates that shocks of equal magnitude have the same influence on the volatility of job loss. Positive shocks likely to create more volatility in the short term than negative shocks if the coefficient is positive. If the coefficient is negative, then negative shocks create more volatility than positive shocks in the short run. Furthermore, when ${\text{isin}}_{\mathrm{t}-1}$ is good or positive news, the overall effect is $\left(1+\gamma \right)/\left|{\text{isin}}_{\mathrm{t}-1}\right|$, but when ${\text{isin}}_{\mathrm{t}-1}$ is bad or negative news, the entire effect is $\left(1-\gamma \right)/\left|{\text{isin}}_{\mathrm{t}-1}\right|$. The EGARCH parameters are not constrained in any manner to ensure that the conditional variance is always positive, whereas the log form of the conditional variance might be negative.

3.3.2. The TGARCH model

The Threshold GARCH (TGARCH) model was proposed by Glosten et al. (1993) and Zokian (1994) and is used to account for leverage effects in financial time series. A TGACRH(1,1) model’s conditional variance specification is as follows:

(10) ${\sigma }_{\mathrm{t}}^2=\omega +\alpha {\text{isin}}_{\mathrm{t}-1}^2+\gamma {\mathrm{d}}_{\mathrm{t}-1}{\text{isin}}_{\mathrm{t}-1}^2+{\beta \sigma }_{\mathrm{t}-1}^2$(10)

where $\omega$, $\alpha$, $\beta$ and $\gamma$ are the volatility long-run averages, prior forecast, symmetric news and negative news, respectively. The indicator function is d, and the dummy variable is ${\mathrm{d}}_{\mathrm{t}-1}$, defined as

(11) ${\mathrm{d}}_{\mathrm{t}-1}=\left\{\begin{array}{c}1\mathrm{i}\mathrm{f}{\text{isin}}_{\mathrm{t}-1}<0,badnews\\ \hspace{0.17em}0\mathrm{i}\mathrm{f}{\text{isin}}_{\mathrm{t}-1}>0,goodnews\end{array}\right.$(11)

The coefficient $\gamma$ represents the asymmetry or leverage term. When $\gamma$ = 0, the TGARCH model departs from the fundamental GARCH model, but when $\gamma \ne 0$, asymmetric shocks exist. When the shock is positive, the effect on volatility is ${\text{isin}}_{\mathrm{t}-1}>0$ (i.e., good news), and when the shock is negative, ${\text{isin}}_{\mathrm{t}-1}<0$ (i.e., bad news), the effect on volatility is $\alpha +\gamma$. Parameter $\left(\alpha +\gamma \right)/\alpha$ is a simple approach to measure the amount of symmetry. As a result, if the shock is both big and positive, negative shocks have a greater influence on volatility than positive shocks. If $\gamma$ were negative, the exact opposite may be true. On the other hand, $\beta$ quantifies clustering in the conditional variance and $\left(\alpha +\beta +\gamma \right)/2$ measures the persistence of shocks on volatility. The parameters $\left(\frac{\alpha +\beta +\gamma }{2}\right)<1$ established a regularity criterion for the presence of the second moment of the TGARCH(1,1) model, according to Ling and McAleer (2002). To provide positive conditional variance, the TGARCH model requires the following parameter constraints: $\omega \ge 0$, $\alpha \ge 0$, $\beta \ge 0$ and $\alpha +\gamma \ge 0$.

3.3.3. The PGARCH model

To deal with asymmetry, Ding et al. (1993) presented the Power GARCH (PGARCH) model. The conditional variance for a PGARCH(1,1) model is given as

(12) ${\sigma }_{\mathrm{t}}^{\mathrm{d}}=\omega +\alpha {\left(\left|{\text{isin}}_{\mathrm{t}-1}\right|+\gamma {\text{isin}}_{\mathrm{t}-1}\right)}^{\mathrm{d}}+{\beta \sigma }_{\mathrm{t}-1}^{\mathrm{d}}$(12)

where d is the power term with $\mathrm{d}>0$ and $\left|\gamma \right|\le 1$. The parameter $\gamma$ is the leverage effect, and when $\mathrm{d}>0,$ $\gamma \ne 0$ are significant, we may say there is asymmetry or a leverage effect. The PGARCH(1,1) model replicates a GARCH(1,1) model for the power term when $\mathrm{d}$ = 2 and $\gamma =0$. If d = 1, the conditional standard deviation will be estimated; and furthermore, the PGARCH model’s effect on volatility is equivalent to that of the TGARCH model.

3.4. Distribution assumption of the error (ε_t)

It is widely accepted that volatile and clustered time series data are not normally distributed. There is excessive kurtosis and a big tail in the residual distribution of the estimated regression. We estimate the GARCH, EGARCH, TGARCH, and PGARCH models in this work by assuming that ${\text{isin}}_{\mathrm{t}}$ has a normal, student-t, and generalised error distribution (G.E.D.). This explains the excess kurtosis and fat tails observed in the time series residuals (see Bollerslev, 1987; Nelson, 1991). These distributions correctly depict the residuals series’ excess kurtosis and skewness. The sensitivity and appropriateness of the assumption findings were seen by changing the distribution assumption from normal to student-t distribution to generalised error distribution.

3.5. Model selection criteria

The best fitting model will be chosen in this study based on the following criteria: (a) diagnostic tests, (b) model selection criteria, and (c) assessment of in-sample and out-of-sample forecasting performances. We use the correlogram Q-statistics (Ljung & Box, 1978) for diagnostic checks and the ARCH LM test (Engle, 1982) for residuals. The Q-statistic (Ljung & Box, 1978) is used to test for serial correlation in the mean equation, whereas the ARCH LM test is used to assess if the variance equation residuals are heteroskedastic. In contrast, we use three model selection criteria to identify the best model: the Akaike information (Akaike, 1974), the Schwarz criterion (Schwarz, 1978), and the Hannan-Quinn criterion (Hannan & Quinn, 1979). Each criterion is based on a probability function, is connected to the others, and may be used separately. The best-fitting model will be the one with the lowest value.

3.6. Forecasting performance measures

In this study, we use three measures to compare the performance accuracy of several competing models: Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), and Theil Inequality Coefficient (Theil, 1967). The model with the smallest forecast error is the most accurate and appropriate.

The following formulas are used to compute the RMSE, MAE, and Theil inequality coefficient,

$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{{\sum }_{\mathrm{t}=1}^{\mathrm{T}}{\left({\sigma }^2-{\hat{\sigma }}^2\right)}^2}{\mathrm{T}}}$, $\mathrm{M}\mathrm{A}\mathrm{E}=\frac{{\sum }_{\mathrm{t}=1}^{\mathrm{T}}\left|{\sigma }^2-{\hat{\sigma }}^2\right|}{\mathrm{T}}$ and $\mathrm{T}\mathrm{h}\mathrm{e}\mathrm{i}\mathrm{l}=\frac{\sqrt{\frac{{\sum }_{\mathrm{t}=1}^{\mathrm{T}}{\left({\sigma }^2-{\hat{\sigma }}^2\right)}^2}{\mathrm{T}}}}{\sqrt{\frac{{\sum }_{\mathrm{t}=1}^{\mathrm{T}}{\sigma }^2}{\mathrm{T}}}+\sqrt{\frac{{\sum }_{\mathrm{t}=1}^{\mathrm{T}}{\hat{\sigma }}^2}{\mathrm{T}}}}$

where T denotes the number of observations, while ${\sigma }^2$ and ${\hat{\sigma }}^2$ denote the actual variance (volatility), and forecasted volatility, respectively. The RMSE computes the difference between the true and estimated values and utilises the sum to assess a model’s prediction potential. According to the criteria, decreasing RMSE values boost the model’s predictive potential. The MAE criteria compute the forecast’s bias and the absolute divergence from the data. While the MAE provides equal weight to deviation over and under estimates, the RMSE gives more weight to substantial forecast errors. The final measure to examine is Theil inequality coefficient, a scale-invariant measure that always lies between zero and one, with zero denoting a perfect fit. This, in turn, occurs precise or without errors.

3.7. Data sources

We utilised real-time administrative data on job loss obtained by the Malaysian Office of Employment Insurance System (EIS) from 1 January 2020 to 31 December 2020. Employees in the formal private sector who lose their jobs are covered under the insured person’s policy when they lose their jobs. On the other hand, we used the total number of COVID-19 new cases and new deaths to represent the coronavirus pandemic. The COVID-19 Government Response Tracker (OxCGRT) database (Hale et al., 2020) was used to gather daily statistics on the number of confirmed new cases and confirmed new deaths, which can be viewed at https://covidtracker.bsg.ox.ac.uk/. All of the variables were converted to logarithms. We employ the formula, $log{\mathrm{x}}_{\mathrm{t}}=\mathrm{l}\mathrm{o}\mathrm{g}\left[{\mathrm{x}}_{\mathrm{t}}+\sqrt{\left({\mathrm{x}}_{\mathrm{t}}^2+1\right)}\right]$ should transform all variables into logarithms (Busse & Hefeker, 2007) in order to keep the sign of ${\mathrm{x}}_{\mathrm{t}}$.

4. Results and discussions

4.1. Data descriptions and preliminary analysis

Table 1 shows the levels and first-difference statistics for the three variables. The level series for job losses has a mean of 6.11, with maximum and lowest values of 8.03 and 2.64, respectively, and a standard deviation of 0.80. The mean, maximum, and lowest values for new cases are 4.30, 8.53, and 0.0, respectively, while the values for new deaths are 0.69, 3.18, and 0.0, respectively. Because the series’ means are positive, we may assume that the series are increasing over time. Both new cases and new deaths had standard deviations of 2.55 and 0.92, respectively. The standard deviation indicates that news about new COVID-19 cases is the most erratic. Job loss, new cases, and log changes in new deaths all have negative skewness, indicating that these series have a longer or fatter tail on the left side of the distribution. It suggests that the series falls more than it increases.

Positive skewness is observed in new deaths, log changes in job loss, and new cases, indicating that these series have a longer or fatter tail on the right side of the distribution and that the series grows more frequently than it declines. On the other hand, all series (with the exception of new cases and deaths in level) have a high kurtosis of larger than 3.0, indicating the existence of fat tails and a leptokurtic series. However, the Jarque-Bera (Jarque & Bera, 1980) and Anderson-Darling (Anderson & Darling, 1952) tests indicate that all variables in the series are non-normal. The null hypothesis of normality may be rejected in all cases at the 1% level of significance, suggesting that none of the three series (both levels and first-differences) follow a normal distribution. The quantile-quantile (Q-Q) plots of each of the series displayed in Figure 4 clearly highlight the non-normality of the series in levels as well as log changes, in agreement with daily data that exhibit skewness and kurtosis in the series.

Table 2 shows the estimated regression findings, which suggest that the relationships are spurious. Because regressing non-stationary variables contradicts the assumption of stationarity of the ordinary least squares estimator, and since the level regression is likely to contain non-stationary variables, thus, resulting in incorrect estimates (Granger & Newbold, 1974). Regressing variables on changes, on the other hand, will result in spurious regression results. The short-run model in this case lacks long-run knowledge on the factors at play. Furthermore, the residuals from both estimated regressions are not normally distributed, indicating that the regression is misspecified with non-constant variance.

Table 2. Regression estimates for loss of employment with COVID-19 new cases and new deaths

Dependent/Independent variables	$\mathrm{j}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{s}{\mathrm{s}}_{\mathrm{t}}$	$\mathrm{\Delta }\mathrm{j}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{s}{\mathrm{s}}_{\mathrm{t}}$
$\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{t}$	6.0145***	0.0025
	(70.168)	(0.0688)
$\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}$	0.0486**
	(2.0663)
$\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}$		−0.0291
		(−0.6508)
$\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}}$	−0.1620**
	(−2.4659)
$\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}}$		−0.0469
		(−0.8325)
${\mathrm{R}}^2$	0.0170	0.0030
SER	0.7919	0.7181
LM${\chi }^2\left(1\right)$	[0.000]	[0.5742]
ARCH${\chi }^2\left(1\right)$	[0.000]	[0.0536]

Notes: Asterisks ***,**,* denote statistically significant at 1%, 5% and 10% level, respectively. LM ${\chi }^2\left(1\right)$ and ARCH${\chi }^2\left(1\right)$ denote the Lagrange multiplier test for serial correlation of order one and heteroscedasticity of order one in the OLS equations, respectively. Figures in round brackets (…) are t-statistics, while figures in square brackets […] are p-values. ${\mathrm{R}}^2$ and SER denote R-squared and standard error of regression, respectively. $\mathrm{j}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{s}{\mathrm{s}}_{\mathrm{t}}$ refers to log in job loss; while $\mathrm{\Delta }\mathrm{j}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{s}{\mathrm{s}}_{\mathrm{t}}$.denotes log changes in job loss. All variables are in logarithm.

According to the classic econometric text book, one of the assumptions in estimating a regression model is that all variables must be stationary. According to research, the bulk of macroeconomic time series variables are stationary in terms of their levels or first-differences (Nelson & Plosser, 1982; Perron, 1988). The first-differences formula for calculating a variable $\mathrm{l}\mathrm{o}\mathrm{g}{\mathrm{x}}_{\mathrm{t}}$ is $\mathrm{\Delta }\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{x}=\mathrm{l}\mathrm{o}\mathrm{g}({\mathrm{x}}_{\mathrm{t}}/{\mathrm{x}}_{\mathrm{t}-1})$. Figure 3 depicts job loss, new cases, and new deaths in levels, as well as log changes (differences) in the job loss, new cases, and new deaths, as well as their unexpected logarithmic oscillations. Figures 3(b), (c), and (d) depict the clustering and volatility in log changes of new cases, deaths, and job loss, respectively. Clustering for COVID-19 new cases is seen in the second and fourth quarters of 2020, whereas clustering for COVID-19 new fatalities occurs in the second and fourth quarters of that year.

Table 2 displays the results of two regressions computed for job loss using new cases and new deaths as regressors. Results in Table 3 clearly suggest that all series are nonstationary variables and that they are I(1). In column 2, of Table 2, we have estimated the regression in level, and in column 3, we report the regression findings in first-differences. Figures 5(a,b) show quantile-quantile (Q-Q) plots of residuals from estimating regressions in level and first-differences, respectively. Regression in log changes, on the other hand, reveals that neither of the two regressors has any influence on job loss. According to the estimated regression findings, both new cases and new deaths are significant. Despite both having relatively low R-squared, regression in changes has a lower standard error of regression. The diagnostic test scores, on the other hand, reveal some surprising findings. While the homoscedastic test failed the regression in changes (at the 10% level), neither the no autocorrelation nor the regression in levels failed. Furthermore, as the Q-Q graphs for both residuals clearly show, the residuals in both regressions are not normally distributed. Our findings clearly suggest the need for more appropriate models, such as Engle’s (1982) Generalized Autoregressive Conditional Heteroscedasticity (GARCH) model and its multiple relatives of GARCH models, which can cope with kurtosis and excess volatility in the variables.

Table 3. Results of unit root tests

Types	Level		First-difference
	Intercept	Intercept+trend	Intercept	Intercept+trend
Job loss	−2.1006 (12)	−3.0452 (12)	−9.0972 (12)***	−9.1143 (12)***
New cases	−1.1899 (4)	−1.7475 (4)	−14.398 (3)***	−14.378 (3)***
New deaths	−2.2756 (4)	−2.9340 (4)	−21.935 (1)***	−21.910 (1)***

Notes: Asterisks *** denotes statistically significant at 1% level. The critical values are referred to MacKinnon (1996); Figures in round bracket (…) denote optimal lag length chosen using SBC criteria.

4.2. Model evaluations

The results of the ARDL(3,3,0)-GARCH(1,1)-M and ARDL(3,3,0)-EGARCH(1,1)-M models for Malaysia’s job losses from January to December 2020 are shown in Table 5. The conditional mean equation, which is displayed in Table 4 as the error-correction model, was derived using the ARDL(3,3,0). Our findings in Table 4 demonstrate unequivocally that heteroscedasticity exists in the residuals of the ECM model, and in light of this, we have introduced the conditional variance, designated by the symbol ${\sigma }_{\mathrm{t}}^2$, as an additional regressor to the conditional mean equation. In order to assess the volatility of job loss in this study, we developed a GARCH-M model.

Table 4. Results of cointegration tests

Dependent/Independent variables	ARDL(3,3,0)	Conditional ARDL-ECM model	Long-run model	Error-correction Model (short-run)
$\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{t}$	2.6648***	2.6648***	5.9215***
	(8.2453)	(8.2453)	(38.472)
$\mathrm{j}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{s}{\mathrm{s}}_{\mathrm{t}-1}$	0.7313***	−0.4500***
	(13.875)	(−8.5025)
$\mathrm{j}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{s}{\mathrm{s}}_{\mathrm{t}-2}$	−0.2851***
	(−4.4952)
$\mathrm{j}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{s}{\mathrm{s}}_{\mathrm{t}-3}$	0.1038**
	(1.9793)
$\mathrm{\Delta }\mathrm{j}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{s}{\mathrm{s}}_{\mathrm{t}-1}$		0.1813***		0.1813***
		(3.4176)		(3.4360)
$\mathrm{\Delta }\mathrm{j}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{s}{\mathrm{s}}_{\mathrm{t}-2}$		−0.1038**		−0.1038**
		(−1.9793)		(−1.9915)
$\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}$	−0.0780*		0.0916**
	(−1.7348)		(2.1109)
$\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-1}$	−0.0752	0.0412**
	(−1.5999)	(2.0817)
$\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-2}$	0.1248***
	(2.6463)
$\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-3}$	0.0698
	(1.5328)
$\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}$		−0.0780*		−0.0780*
		(−1.7348)		(−1.7527)
$\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-1}$		−0.1946***		−0.1946***
		(−3.7858)		(−3.8903)
$\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-2}$		−0.0698		−0.0698
		(−1.5328)		(−1.5474)
$\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}}$	−0.1211**		−0.2692**
	(−2.2455)		(−2.2841)
$\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{d}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{h}{\mathrm{s}}_{\mathrm{t}-1}$		−0.1211**
		(−2.2455)
$\mathrm{E}\mathrm{C}{\mathrm{M}}_{\mathrm{t}-1}$				−0.4500***
				(−8.5643)
${\mathrm{R}}^2$	0.4275	0.2847	-	0.2883
LM${\chi }^2\left(1\right)$	[0.5692]	-	-	[0.7519]
ARCH${\chi }^2\left(1\right)$	[0.0335]	-	-	[0.0248]
Bound F-statistic	-	18.182***	-	-

Notes: Asterisks ***,**,* denote statistically significant at 1%, 5% and 10% level, respectively. Figures in round brackets (…) are t-statistics, while figures in square brackets […] are p-values. ${\mathrm{R}}^2$ denotes R-squared. $\mathrm{j}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{s}{\mathrm{s}}_{\mathrm{t}}$ refers to job loss. LM ${\chi }^2\left(1\right)$ and ARCH${\chi }^2\left(1\right)$ denote the Lagrange multiplier test for serial correlation of order one and heteroscedasticity of order one in the ARDL equation, respectively ARDL(p,q1,q2) denotes optimal lag length chosen using AIC criteria.

Table 5. Parameter estimates of ARDL(3,3,0)-GARCH(1,1)-M and ARDL(3,3,0)-EGARCH(1,1)-M

Independent variables	ARDL(3,3,0)-GARCH(1,1)-M:			ARDL(3,3,0)-EGARCH(1,1)-M:
	Normal	Student’s t	GED	Normal	Student’s t	GED
Mean equation
${\delta }_0\hspace{0.17em}\left(\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{t}\right)$	1.0537***	1.0472	0.8440	−0.2274***	0.7190***	0.2130***
	(2.5872)	(1.5831)	(1.6096)	(−3.4730)	(3.1519)	(3.0880)
$\pi \mathrm{E}\mathrm{C}{\mathrm{M}}_{\mathrm{t}-1}$	−0.3845***	−0.3608***	−0.2974***	−0.4035***	−0.0543**	−0.0565**
	(−5.5116)	(−6.6739)	(−7.8204)	(−8.7297)	(−2.1631)	(−2.5519)
${\delta }_1\mathrm{\Delta }\mathrm{j}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{s}{\mathrm{s}}_{\mathrm{t}-1}$	0.1759***	0.1153**	0.0360	−0.0934**	0.1619*	0.0900**
	(2.3363)	(2.0842)	(0.9431)	(−2.4910)	(1.7642)	(1.9750)
${\delta }_2\mathrm{\Delta }\mathrm{j}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{s}{\mathrm{s}}_{\mathrm{t}-2}$	−0.1364**	−0.1029**	−0.0062	−0.1001***	−0.1211**	0.0285
	(−2.2593)	(−2.1387)	(−0.2204)	(−3.8742)	(−2.5704)	(0.9398)
${\delta }_3\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}$	−0.0812**	−0.0680*	−0.1016***	−0.0836***	−0.0948***	−0.0937***
	(−2.0340)	(−1.8561)	(−3.6049)	(−3.5104)	(−2.7601)	(−4.8875)
${\delta }_4\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-1}$	−0.1770***	−0.1557***	−0.1595***	−0.1629***	−0.1504***	−0.1088***
	(−3.7755)	(−3.8670)	(−5.0559)	(−7.8241)	(−4.3637)	(−4.8202)
${\delta }_5\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-2}$	−0.0651	−0.0820**	−0.0682**	−0.0913***	−0.0101	−0.0301
	(−1.2857)	(−1.9703)	(−2.4835)	(−4.5609)	(−0.2516)	(−1.3809)
${\lambda \mathrm{l}\mathrm{o}\mathrm{g}\sigma }_{\mathrm{t}}^2$	1.0035***	1.0438	0.8882	−0.1650***	0.5994***	0.2179***
	(2.6340)	(1.5250)	(1.6071)	(−4.2421)	(3.0789)	(3.7818)
Variance equation
$\mathrm{c}$ (constant)	0.1884	0.1889**	0.2370***	−0.3508***	−0.7744***	−0.4843***
	(1.5693)	(2.0965)	(3.4193)	(−3.5821)	(−3.9977)	(−3.6169)
$\alpha$ (ARCH effect)	0.0542*	0.0542	0.0658	−0.6828***	0.1026*	0.0301
	(1.6501)	(1.2259)	(1.3543)	(−5.7886)	(1.7633)	(1.2895)
$\beta$ (GARCH effect)	0.4153	0.4639**	0.3384**	0.3026***	0.4032***	0.5813***
	(1.1694)	(2.0538)	(2.2615)	(4.0766)	(3.0049)	(5.9841)
$\gamma$ (Leverage effect)				−0.7794***	−0.4119***	−0.5990***
				(−9.3706)	(−4.2888)	(−5.2422)
$\alpha +\beta$	0.4695	0.5181	0.4042	−0.3802	0.5058	0.6114
${\mathrm{R}}^2$	0.2988	0.2893	0.2562	0.2702	0.2850	0.1689
SER	0.6053	0.6094	0.6234	0.6175	0.6112	0.6590
Q-stat:	[0.9480]	[0.3880]	[0.5690]	[0.0940]	[0.1810]	[0.4230]
ARCH test:	[0.6456]	[0.6914]	[0.2761]	[0.9136]	[0.1299]	[0.0164]
Criteria:
AIC	1.8607	1.7768	1.7670	1.6066	1.6717	1.6412
SC	1.9792	1.9061	1.8962	1.7359	1.8118	1.7812
HQC	1.9078	1.8282	1.8184	1.6580	1.7274	1.6968

Notes: Notes: Asterisks ***,**,* denote statistically significant at 1%, 5% and 10% level, respectively. Figures in round brackets (…) are z-statistics, while figures in square brackets […] are p-values. ${\mathrm{R}}^2$ and SER denote R-squared and standard error of regression, respectively; while Q-stat and ARCH test are test for autocorrelation and heteroscedasticity, respectively. AIC, SC and HQ denote Akaike information criteria, Schwarz criteria and Hannan-Quinn criteria, respectively.

A similar technique was used by Holmes and Maghrebi (2016) and Haughton and Iglesias (2017) to provide an ECM model for the conditional mean equation. However, Holmes and Maghrebi’s (2016) analysis on the link between the US stock market and the unemployment rate failed to uncover any cointegration relationship, forcing them to estimate a short-run mean equation without the $\mathrm{E}\mathrm{C}{\mathrm{M}}_{\mathrm{t}-1}$ component. Haughton and Iglesias (2017) utilised the ARDL-GARCH(1,1) model to investigate the links between exchange and the stock market in the Caribbean and Latin America. They described the conditional mean equation in the mean equation using the unconstrained error-correction model.

Table 5 includes the GARCH-M (columns 2–4) and EGARCH-M (columns 5–7) models, as well as three types of error distributions: normal, student-t, and generalised. The table displays the mean and variance equations, goodness of fit (R-squared), standard error of regression (SER), Q-statistics for the serial correlation test, ARCH test for the presence of heteroscedatic in the error, and three model section criteria—AIC, SC, and HQC. Our GARCH-M model’s mean equation suggests that virtually all measured parameters are important. At the 1% level, the calculated coefficients for the error-correction component $\pi$, are negative and significant. The relevance of $\pi$ shows a cointegration or long-term link between job loss and new COVID-19 cases and deaths.

The short-run variables—changes in lagged (one and two period) job loss, and current and lagged (one and two period) new cases of COVID-19 all had an effect on the current changes in job loss. When the error distribution is normal, the conditional variance or volatility in the job loss affects only the most recent changes in the job loss. In the variance equation, however, the GARCH impact is positive and significant in the student-t and generalised error distribution models, but the ARCH effect is only positive and significant in the normal error distribution model. The sum of the ARCH and GARCH effects, which varies from 0.40 to 0.52 and is less than one, suggests that the volatility in the job loss is relatively moderately persistence. Furthermore, the ARDL(3,3,0)-GARCH(1,1)-M models lack serial correction and heteroscedasticity.

The estimated EGARCH-M model, on the other hand, shows that the error-correction model, in which the parameter of the $\mathrm{E}\mathrm{C}{\mathrm{M}}_{\mathrm{t}-1}$ component is negative and significant at the 1% level, fits the data quite well. Changes in the present, lagged one and two periods in new cases of Covid-19, as well as changes in the lagged one and two periods of job loss, are significant short-run variables with negative signs in the majority of cases. However, in the normal error distribution model, volatility in job loss is large and has a negative sign, but it has a positive influence on fluctuations in job loss in both the student-t- and generalised error distribution models. For the normal error distribution, positive volatility in the job loss implies that higher volatility will increase current changes in the job loss, while negative volatility in the job loss implies that higher volatility will decrease current changes in the job loss.

With the exception of the constant component, the ARCH and GARCH influences on the variance equation in the student-t error distribution model are both positive and significant. The total of the ARCH and GARCH effects is less than one for both the student-t and generalised error distribution models, and it equals 0.5 for the student-t error distribution model, indicating that volatility persistency is mild. The estimated parameter, which is negative and significant at the 1% level, demonstrates the influence of leverage (or asymmetry) on changes in job loss. Bad news shocks, according to the negative leverage effect, create greater volatility than good news shocks (good news). With the exception of the model for the generalised error distribution, the diagnostic tests show that the ARDL(3,3,0)-EGARCH(1,1)-M model is devoid of serial correlation and heteroscedasticity.

The TGARCH-M and PGARCH-M models, like the EGARCH-M model, can test for asymmetry in the volatility of job loss. Table 6 shows the results of estimating both of these models. We evaluated each model using one of three error distributions: the normal, student-t, and generalised error distributions. TGARCH-M and PGARCH-M estimated conditional mean equations show that the data match the ECM estimated equations quite well. The $\mathrm{E}\mathrm{C}{\mathrm{M}}_{\mathrm{t}-1}$ term is negative and significant at the 1% level, indicating cointegration between job loss and new cases and deaths. In the short-run, the lagged one and two periods in job loss, as well as the present and lagged one and two periods in new cases Covid-19 influence the current changes in job loss. The mean equation also revealed that the volatility of job loss influences the present changes in job loss.

Table 6. Parameter estimates of ARDL(3,3,0)-TGARCH(1,1)-M and ARDL-PGARCH(1,1)-M

Independent variables	ARDL(3,3,0)-TGARCH(1,1)-M:			ARDL(3,3,0)-PGARCH(1,1)-M:
	Normal	Student’s t	GED	Normal	Student’s t	GED
Mean equation
${\delta }_0\hspace{0.17em}\left(\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{t}\right)$	0.0803	0.4747**	0.6918***	−0.1023**	0.6587	0.6140***
	(0.5306)	(2.4329)	(2.7632)	(−2.2053)	(1.4117)	(4.3756)
$\pi \mathrm{E}\mathrm{C}{\mathrm{M}}_{\mathrm{t}-1}$	−0.2492***	−0.1576***	−0.1143***	−0.3901***	−0.4725***	−0.0701**
	(−4.0611)	(−3.5593)	(−3.8741)	(−9.2235)	(−6.1104)	(−2.3634)
${\delta }_1\mathrm{\Delta }\mathrm{j}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{s}{\mathrm{s}}_{\mathrm{t}-1}$	0.2549***	0.1811**	0.0828	0.0700***	0.1639***	−0.0253
	(3.8225)	(2.4954)	(1.4921)	(4.2602)	(3.0377)	(−0.5276)
${\delta }_2\mathrm{\Delta }\mathrm{j}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{s}{\mathrm{s}}_{\mathrm{t}-2}$	−0.1511***	−0.0995**	−0.0010	−0.1118***	−0.0597	−0.0374
	(−2.9828)	(−2.1611)	(−0.0414)	(−3.8089)	(−1.1557)	(−1.0865)
${\delta }_3\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}}$	−0.0641*	−0.0564*	−0.0980***	−0.0383	−0.0645*	−0.0701***
	(−1.7951)	(−1.7217)	(−4.3801)	(−1.6435)	(−1.7730)	(−3.0091)
${\delta }_4\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-1}$	−0.1308***	−0.1128***	−0.1245***	−0.1294***	−0.1605***	−0.1030***
	(−3.2193)	(−3.1794)	(−5.3889)	(−7.7116)	(−3.9082)	(−3.7481)
${\delta }_5\mathrm{\Delta }\mathrm{n}\mathrm{e}\mathrm{w}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}{\mathrm{s}}_{\mathrm{t}-2}$	−0.0343	−0.0177	−0.0112	−0.0372	−0.0734*	−0.0102
	(−0.7491)	(−0.4287)	(−0.4609)	(−1.2953)	(−1.8667)	(−0.3801)
${\lambda \mathrm{l}\mathrm{o}\mathrm{g}\sigma }_{\mathrm{t}}^2$	0.0811	0.4150***	0.6740***	−0.0932**	0.6917	0.4905***
	(0.6723)	(2.5804)	(2.9043)	(−2.2564)	(1.4190)	(4.4770)
Variance equation
$\mathrm{c}$ (constant)	0.1548***	0.1962***	0.2038***	0.5394***	0.0984***	0.3666***
	(10.364)	(3.4248)	(4.3116)	(6.8302)	(10.169)	(8.8817)
$\alpha$ (ARCH effect)	−0.1452***	−0.1005**	−0.0333	0.2081***	0.0299*	0.1996***
	(−10.669)	(−2.3643)	(−1.1723)	(3.7644)	(1.7022)	(3.4812)
$\beta$ (GARCH effect)	0.3899***	0.2883**	0.3453***	0.2301**	0.8529***	0.3350***
	(41.094)	(2.2572)	(3.4170)	(2.3569)	(138.43)	(3.8049)
$\gamma$ (Leverage effect)	0.6966***	0.6257**	0.3066**	0.9372***	−0.3995	0.9691***
	(4.1254)	(2.4138)	(2.1725)	(12.446)	(−0.7352)	(10.677)
$\mathrm{d}$ (Power)				0.2343*	0.5462	0.6948***
				(1.6496)	(1.4317)	(3.8487)
$\alpha +\beta$	0.2447	0.1878	0.3120	0.4382	0.8828	0.5346
${\mathrm{R}}^2$	0.2563	0.2985	0.2583	0.2739	0.2832	0.2505
SER	0.6234	0.6054	0.6225	0.6160	0.6120	0.6258
Q-stat:	[0.1770]	[0.5430]	[0.4100]	[0.1110]	[0.7940]	[0.0430]
ARCH test:	[0.1726]	[0.1021]	[0.0634]	[0.0525]	[0.1114]	[0.1892]
Criteria:
AIC	1.7763	1.7176	1.6942	1.7395	1.7987	1.6725
SC	1.9056	1.8576	1.8343	1.8796	1.9495	1.8233
HQC	1.8277	1.7732	1.7499	1.7952	1.8586	1.7324

Notes: Notes: Asterisks ***,**,* denote statistically significant at 1%, 5% and 10% level, respectively. Figures in round brackets (…) are z-statistics, while figures in square brackets […] are p-values. ${\mathrm{R}}^2$ and SER denote R-squared and standard error of regression, respectively; while Q-stat and ARCH test are test for autocorrelation and heteroscedasticity, respectively. AIC, SC and HQ denote Akaike information criteria, Schwarz criteria and Hannan-Quinn criteria, respectively.

The TGARCH-M model results indicate that the ARCH effects are negative while the GARCH effects are positive; nevertheless, the leverage impact, which is positive and significant, clearly indicates the presence of an asymmetric effect of volatility on job loss. The importance of the leverage effect implies that negative shocks (bad news) have a higher influence on conditional variance (volatility) than positive shocks (good news) of equal size. In addition, the TGARCH-M model cannot reject the null hypothesis of non-serial correlation in the mean equation and homoscedastic residuals in the variance equation.

The PGARCH-M model, on the other hand, shows that both the ARCH and GARCH effects are positive and significant. The combined ARCH and GARCH impacts are less than one. In the models of the normal and generalised error distributions, the leverage effect is strong and positive. Positive shocks, according to the positive leverage, are associated with more volatility than negative shocks. The power term d is positive and significant in the models for the normal and generalised error distributions. Because $\mathrm{d}\ne 2$, we may conclude that the model is not a standard GARCH model. Contrary to the Q-statistic test for serial correlation, which indicates that only model estimates using the normal and Student’s t-error distributions are free of serial correlation, the ARCH test for heteroscedasticity indicates that none of the three estimated variance equations exhibit heteroscedastic error.

4.3. Model forecast accuracy

We used four different GARCH-M models to explain for the variability in job loss in Malaysia during the COVID-19 pandemic from January to December 2020. The next stage is to determine which of the four models best depicts the volatility in job loss, each with three distinct error distribution assumptions. When choosing the best model for this study, the following criteria will be considered: (1) model selection criteria; (2) estimated parameters that suit the theory; (3) in-sample forecasting capability; and (4) out-of-sample forecasting ability. According to the results in Tables 5 and 6, the error distribution models with the lowest AIC, SC, and HQC are the generalised error distributions for GARCH-M, TGARCH-M, and PGARCH-M, as well as the normal error distribution for EGARCH-M.

However, only the GARCH effect is significant for the GARCH-M model, but the non-negativity criteria of the EGARCH-M model was not fulfilled since the ARCH effect has a negative sign and the sum of the ARCH and GRACH effects is similarly negative. In contrast, the G.E.D model for the TGARCH-M model failed the non-negativity condition when the ARCH effect was negative, but the G.E.D model for the PGARCH-M model fulfilled all requirements but revealed heteroscedasticity in the variance equation residuals.

Nonetheless, Table 7 shows how well each model predictions with all three alternative residual distribution assumptions. Regardless of the model chosen based on the aforementioned model selection criteria, the optimum model is G.E.D for GARCH-M; student-t for EGARCH-M; and normal error distribution for both TGARCH-M and PGARCH-M models. However, the normal error distribution for the GARCH-M and TGARCH-M models, as well as the G.E.D for the EGARCH-M and PGARCH-M models, have the least RMSE, MAE, and Theil inequality coefficients. However, because both the GARCH-M and TGARCH-M models failed to fulfil the non-negativity requirements, the EGARCH-M and PGARCH-M models are the only ones that may be used.

Table 7. Evaluation of in-sample and out-of-sample forecasting for (log) loss of employment

Models	In-sample forecasting:			Out-of-sample forecasting:
	RMSE	MAE	Theil	RMSE	MAE	Theil
ARDL(3,3,0)-GARCH(1,1)-M:
Normal	11.077	9.5431	0.4857	1.5093	1.1759	0.1214
Student-t	21.050	20.222	0.6346	2.2873	1.6472	0.1708
G.E.D.♠	9.8860	9.0964	0.4520	1.7008	1.2738	0.1334
ARDL(3,3,0)-EGARCH(1,1)-M:
Normal♠	7.9783	7.2856	0.4839	1.4866	1.2614	0.1248
Student-t	4.8027	4.5479	0.2834	0.8902	0.6964	0.0716
G.E.D.	5.8756	4.7158	0.5713	0.6621	0.5371	0.0551
ARDL(3,3,0)-TGARCH(1,1)-M:
Normal	4.8022	3.5811	0.3064	1.0855	0.8670	0.0890
Student-t	19.200	17.676	0.6164	1.8321	1.4039	0.1388
G.E.D. ♠	10.421	9.9366	0.4623	1.1710	0.9052	0.0939
ARDL(3,3,0)-PGARCH(1,1)-M:
Normal	7.5315	5.6032	0.4192	1.8989	1.7559	0.1666
Student-t	22.223	21.191	0.6480	2.3617	1.6665	0.1752
G.E.D. ♠	12.288	11.314	0.5064	1.3290	1.0608	0.1082

Notes: RMSE, MAE and Theil refer to root mean square error, mean absolute error, mean absolute percent error and Theil inequality coefficient, respectively.

How can we choose the best model from this complex analysis? We do this by using an elimination procedure. Because the ARCH effect does not fulfil the non-negativity condition of the EGARCH-M model, for example, the normal error distribution model can be ignored for the EGARCH-M model, but the G.E.D model failed to satisfy the variance equation’s homoscedastic error property. As a result, the student-t error distribution is the best model for the EGARCH-M model. In fact, when compared to the G.E.D model, this model has the second-lowest RMSE, MAE, and Theil inequality coefficients, as well as the lowest total RMSE, MAE, and Theil inequality coefficients.

In-sample forecasting, on the other hand, may be used to determine if the PGARCH-M model is the best. In this regard, it is obvious that, despite having the lowest RMSE, MAE, and Theil inequality coefficients for out-of-sample forecasting accuracy, the G.E.D. model displays heteroscedastic error in the variance equation. In terms of predicting abilities, the normal error distribution model has the second-lowest RMSE, MAE, and Theil inequality coefficients and the smallest RMSE, MAE, and Theil inequality coefficients overall. Based on the information supplied, the normal error distribution model is the best fit for PGARCH-M. However, after comparing EGARCH-M (Student-t error distribution) and PGARCH-M (normal error distribution), it is plausible to infer that the ARDL(3,3,0)-EGARCH(1,1)-M is the best model to describe volatility in the job loss in Malaysia from January to December 2020. The RMSE, MAE, and Theil inequality coefficients reveal that the EGARCH-M model has smaller in-sample and out-of-sample forecasting accuracy than the PGARCH-M model. These findings are consistent with recent studies (Adhikari et al., 2021; Arndt et al., 2020; Asfaw, 2021; Saari et al., 2022; Shafi et al., 2020).

5. Conclusion

The labour market is one of the economic areas most affected by the COVID-19 pandemic. According to surveys, the unemployment rate soared in a number of nations shortly after the shutdown. Similarly, as the Malaysian government took lockdown measures to combat the COVID-19 outbreak, the number of individuals losing their jobs increased. In this study, we investigated the long-run and short-run effects of new COVID-19 cases and deaths on the Malaysian labour market. Using daily data on job loss, the number of new cases, and deaths, we found that job loss had long-run associations with the number of new cases and the number of new deaths. The COVID-19 pandemic measures will have a long-term impact on Malaysian job losses from January to December 2020.

In the near run, we detect unpredictability and clustering in job losses during the COVID-19 pandemic in Malaysia. To simulate the volatility and clustering in changes in job loss, we utilise the symmetric GARCH-M model and three asymmetric GARCH-M model variations, namely the EGARCH-M, TGARCH-M, and PGARCH-M models. The asymmetric GARCH models will capture the leverage effects, implying that the volatility in the job loss series will be more volatile in reaction to negative shocks (bad news) than in response to positive shocks (good news). Based on model selection criteria and in-sample and out-of-sample forecasting performances, the current study reveals that the EGARCH-M model is the best model to explain the volatility and clustering in job loss in Malaysia over the 2020 pandemic period. In Malaysia, the variance in job loss is impacted by short- and long-term unemployment statistics, the number of new cases, the number of new deaths, and the conditional variance. As a result, a short-run model that considers both short- and long-term information can provide a better model and be more relevant for projecting job losses in Malaysia.

The implications of this study are as follows: (i) the surge in the loss of employment associated with lockdown measures suggest that continues imposition of the lockdown measures will not be long-run solution. Policy intervention such as providing reliable vaccines enable the gradual opening up of the economy and reduce vulnerability arising from lockdown induced job losses. (ii) Since job loss is more volatile in reaction to negative news such as increase in the number of new cases and new deaths, the policy implication is that there is a need for continuous awareness for individual preventive measures rather than lockdowns and the need to caution media houses from over-exaggerating outbreak information.

Our analysis has a practical purpose. One important application is the “value-at-risk” or VaR analysis. Any government or private institutions would want to measure the risks faced in producing outputs. For example, a 1% value at risk is defined as the number of Ringgit that one can be 99% certain exceeds any output losses for the next day. The VaR analysis has been frequently applied by many banks and other financial institutions to measure losses due to the volatility in their portfolios (Engle, 2001).

Lastly, how robust is our results is a question to be answered. However, we leave this issue for the next agenda. Robustness can be tested by using different samples of study, by adding new variables, and also, by testing for outliers. The work by Akpan et al. (2021) has employed the Outlier Modified GARCH-type models. Akpan et al. (2021) noted that “The results obtained revealed that the Outlier Modified GARCH-type models chosen were adequate and sufficiently reducing the value of excess kurtosis in close proximity to the theoretical value. Therefore, the modification of existing GARCH-type models by subtracting the effects of outliers seems to show a substantive improvement in the portion of excess kurtosis captured and thus proves that the Outlier Modified GARCH-type models make difference with the existing ones.”

Disclosure statement

No potential conflict of interest was reported by the author(s).