In this paper, we first derive the energy consumption function equation from the general firm cost function equation, which draws from Fisher-Vanden et al. (2004) and Yao and Yu (2012). The decision model of energy CO₂ emission intensity is then derived. In general, we can write the cost function of the firm as follows: C = c Q , A , P K , P L , P E , (1) where C is the cost input factor, P is the price input factor, K is the capital input factor, L is the labor input factor, E is the energy input factor, Q is the gross output, and A is the production technique level. Assuming that the production function is Cobb–Douglas production function with the Hicks-neutral effect, we can write the cost production function as follows: C = α K + α L + α E α K α K α L α L α E α E − 1 α K + α L + α E A − 1 Q P K α K P L α L P E α E 1 α K + α L + α E , (2) where α _i is the output elasticity of factor i. According to Shephard’s lemma, the demand for a given factor of production is the partial derivative of the cost of the production function concerning the price of such a production factor. Thus, for a determined total output and other factors of production, the equation for the energy consumption function is E = α E A − 1 P K α K P L α L P E α E Q / P E . (3)

Equation 3 shows that the total output value, production technology level, energy price, and total factor input price jointly determine the energy consumption in actual production.

Assumption: P Q = P K α K P L α L P E α E , (4) where α K + α L + α E = 1 . Substituting Eq. 4 into Eq. 3 yields a more intuitive equation for the energy consumption function: E = α E A − 1 Q R P E / Q = α E A − 1 R P , (5) where R P = P K α K P L α L P E α E P E − 1 denotes the price of energy relative to the price of all factors.

Since fossil energy consumption is the primary source of CO2 emissions, that is, CO2 emissions are related to energy consumption and energy structure, we further obtain the determining equation of CO2 emission intensity according to Eq. 5. C I = α E A − 1 R P E S , (6) where CI is the CO₂ emission intensity and E.S. is the energy consumption structure.

Therefore, Eq. 6 shows that under constant returns to scale, the output elasticity of energy and other factors, production technology level, energy structure, and the relative price of energy determine the CO2 emission intensity of enterprises.

As mentioned earlier, this paper mainly examines the effect of green technology innovation on CO₂ emission intensity. At the same time, in the analysis of factors influencing macro CO₂ emission intensity, in addition to scientific and technological innovation, industrial structure and the foreign direct investment have a non-negligible effect on total factor productivity and CO₂ emission. Therefore, we interpret the production technology level in Eq. 6 as the generalized productivity, including scientific and technological innovation, industrial structure, and the foreign direct investment. Among them, to distinguish the difference between the effect of green invention patents GINV and green utility model patents GUTY on CO₂ emission intensity, both of them are included in the analysis framework at the same time, so the formula is A = A 0 exp ⁡   δ l n G I N V + τ l n G U T Y + η l n I S + λ l n F D I , (7) where A ₀ represents the initial value of the technology level, which is a constant. GINV is the green technology innovation level expressed by green invention patents. GUTY is the green technology innovation level expressed by green utility model patents. I.S. represents the industrial structure, and FDI represents the foreign direct investment. δ, τ, η, and λ is the impact coefficients, respectively. Combining Eqs 6, 7, the benchmark model used in the analysis of this paper can be determined as follows: l n C I t = α 0 + α 1 l n G I N V t + α 2 l n G U T Y t + α 3 l n I S t + α 4 l n F D I t + α 5 l n R P t + α 6 l n E S t . (8)

In order to test the long-term cointegration and short-term dynamic relationship between green invention patent GINV, green utility model patent GUTY, industrial structure I.S., foreign direct investment FDI, relative energy price R.P., energy consumption structure E.S., and CO2 emission intensity CI, this paper uses the ARDL model proposed by Charen and Deadman (1997) and developed by Pesaran and Pesaran, 1997; Pesaran and Shin, 1999, to analyze. The method does not require the variables to be homogeneous and single integers, and the analysis model can contain a single integer I (0) of order 0 or I (1) of order 1. This method can still obtain unbiased estimates for small samples or endogenous variables. Combined with the benchmark model Eq. 8, the ARDL (m0, m1, m2..., m6) model is as follows: Δ ln ⁡ C I t = ψ ω   t + ∑ j = 0 m 0 β 0 j Δ ln ⁡ C I t − j + ∑ j = 0 m 1 β 1 j Δ ln ⁡ G I N V t − j + ∑ j = 0 m 2 β 2 j Δ ln ⁡ G U T Y t − j + ∑ j = 0 m 3 β 3 j Δ ln ⁡ I S t − j + ∑ j = 0 m 4 β 4 j Δ ln ⁡ F D I t − j + ∑ j = 0 m 5 β 5 j Δ ln ⁡ R P t − j + ∑ j = 0 m 6 β 6 j Δ ln ⁡ E S t − j + μ t . (9)

In the aforementioned equation, ω _t is a deterministic variable including constant, time, and time trend. Furthermore, the lag order (m ₀ , m ₁ , m ₂ , ..., m ₆ ) is obtained based on AIC, SIC, H.Q., and Adj-R².

Based on Eq. 9, the conditional error correction model (conditional error correction model) is as follows: Δ ln ⁡ C I t = ϕ ω   t + ∑ j = 1 m 0 α 0 j Δ ln ⁡ C I t − j + ∑ j = 0 m 1 α 1 j Δ ln ⁡ G I N V t − j + ∑ j = 0 m 2 α 2 j Δ ln ⁡ G U T Y t − j + ∑ j = 0 m 3 α 3 j Δ ln ⁡ I S t − j + ∑ j = 0 m 4 α 4 j Δ ln ⁡ F D I t − j + ∑ j = 0 m 5 α 5 j Δ ln ⁡ R P t − j + ∑ j = 0 m 6 α 6 j Δ ln ⁡ E S t − j + δ 0 ⁡ ln ⁡ C I t − 1 + δ 1 ⁡ ln ⁡ G I N V T t − 1 + δ 2 ⁡ ln ⁡ G U T Y t − 1 + δ 3 ⁡ ln ⁡ I S t − 1 + δ 4 ⁡ ln ⁡ F D I t − 1 + δ 5 ⁡ ln ⁡ R P t − 1 + δ 6 ⁡ ln ⁡ E S t − 1 + μ t . (10)

Based on Eq. 10, two methods, the Wald test and the t-test, can be applied to examine the cointegration relationship between variables. Among them the Wald test is suitable to examine the joint significance level of each lagged variable with the original hypothesis H 0 : δ 0 = δ 1 = δ 2 = ⋯ = δ n = 0 , while the t-test is applied to examine the significance of the lagged terms of the explained variables with the original hypothesis H 0 : δ 0 = 0 . Both indicate that there is no cointegration relationship between the variables. Pesaran et al. (2001) gave two sets of critical values for all variables obeying I (0) or I (1), respectively: lower and upper critical values. In the specific test, if the value of the F-statistic (or t-statistic) obtained using Eq. 10 is greater (or less) than the upper critical value, it indicates that there is a cointegration relationship between the variables; if it is less (or greater) than the lower critical value, then there is no cointegration relationship.

If the test results show a cointegration relationship, the corresponding long-run elasticity coefficients for each influencing factor can be available via Eq. 10: ϕ ′ = ϕ / 1 − ∑ j = 1 m 0 β 0 j ,   ρ i = ∑ j = 0 m 0 β i j / 1 − ∑ j = 1 m 0 β 0 j   i = 1,2 , . . . , 6 . (11)

This equation leads to the long-run cointegration equation between CO₂ emission intensity and each explanatory variable: e c m   t = ln ⁡ C I t − ϕ ′ ω t + ρ 1 ⁡ ln ⁡ G I N V T   t + ρ 2 ⁡ ln ⁡ G U T Y   t + ρ 3 ⁡ ln ⁡ I S   t + ρ 4 ⁡ ln ⁡ F D I t + ρ 5 ⁡ ln ⁡ R P   t + ρ 6 ⁡ ln ⁡ E S   t . (12)

Furthermore, a linear transformation of Eq. 9 yields the ARDL-ECM model for short-term effects estimation. Δ ln ⁡ C I t = ϕ ω   t + ∑ j = 1 m 0 − 1 λ 0 j Δ ln ⁡ C I t − j + ∑ j = 0 m 1 − 1 λ 1 j Δ ln ⁡ G I N V t − j + ∑ j = 0 m 2 − 1 λ 2 j Δ ln ⁡ G U T Y t − j + ∑ j = 0 m 3 − 1 λ 3 j Δ ln ⁡ I S t − j + ∑ j = 0 m 4 − 1 λ 4 j Δ ln ⁡ F D I t − j + ∑ j = 0 m 5 − 1 λ 5 j Δ ln ⁡ R P t − j + ∑ j = 0 m 6 − 1 λ 6 j Δ ln ⁡ E S t − j − θ e c m t − 1 + μ t . (13)

The coefficient θ error correction term ecm _t−1 reflects the strength of the adjustment of CO₂ emission intensity to its deviation from the long-term equilibrium.

Although the ARDL method does not require the variables to be of the same order, we can analyze either the zeroth order I(0) or the first order I(1). However, for series I(2) and above, the F-statistic and t-statistic of the boundary cointegration test may be deviated (Pesaran et al., 2001). Therefore, carrying out a unit root test for each variable is necessary before using the ARDL method for the cointegration test. Since the sample number of this paper is only 28 years, considering that the DF-GLS (GLS transformed Dickey–Fuller) unit root test method can overcome the defects of low efficacy of commonly used ADF(augmented Dickey–Fuller) unit root test method and P.P. (Phillips–Perron) unit root test method in the small sample unit root test. Therefore, the DF-GLS unit root test method is adopted, and the null hypothesis is that there is a unit root in the sequence. The results in Table 1 show that the lag order selection is based on SIC criterion (Schwarz Information Criterion) and AIC criterion (Akaike Information Criterion), and the maximum lag order is 5.

As can be seen from Table 1, all variables adopted in this paper are first-order integrated I(1) processes, so the ARDL model can be further used to analyze whether there is a long-term cointegration relationship among variables.

Before using Eq. 10 for the cointegration test, we should first calculate the optimal lag order of Eq. 9. We choose the maximum final order as 2, through the empirical sample size of this paper, and then according to AIC, SIC, and H.Q. criteria, the optimal ARDL models without and with trend terms are ARDL (2,2,2,2) and ARDL (2,2,2,1,1,1,2), respectively. Table 2 shows the corresponding estimation results.

According to Table 2, the trend term of ARDL with trend term (2,2,2,1,1,1,2) is not significant. Meanwhile, based on statistical indicators such as AIC and H.Q., it is found that ARDL without a trend term (2,2,2,2,2) is better than the model with a trend term. Therefore, we analyzed the ARDL (2,2,2,2,2) model without trend term as the basis for analysis.

Table 3 shows the results of the ARDL (2,2,2,2,2) cointegration test without a trend term. Considering that the critical value of F-statistic provided by Pesaran et al. (2001) is calculated based on large samples, and the number of samples in this paper is only 28, the critical value obtained by random simulation of six explanatory variables and 30 samples in Narayan and Smyth, 2005 is adopted for F-statistic. The results show that the null hypothesis is rejected at a 1% significance level for both F- and t-statistics. Therefore, in the sample period, there is a long-term stable relationship between CO₂ emission intensity (CI) and green invention patents (GINV), green utility model patents (GUTY), industrial structure (I.S.), foreign direct investment (FDI), relative energy price (R.P.), and energy consumption structure (E.S.).

Based on the boundary cointegration test results, the long-term relationship coefficients between variables are estimated using Eq. 11 with electricity dependence as the explanatory variable. We estimate the short-term dynamic coefficients between variables using Eq. 13; Table 4 shows the results of the long-term and short-term coefficient estimates and diagnostic tests. The results show that the model residuals satisfy the normal distribution and eliminate the phenomena of serial autocorrelation and heteroskedasticity. In contrast, the results of the Ramsey RESET test for the model form misspecification show that the model is correctly set.

Table 4 shows that the error correction term ecm (−1) is significant at the 1% significance level with a negative sign. The larger the absolute value of this coefficient, the faster the system returns to equilibrium after a shock. The regression result of the error correction term is −0.8831, which indicates that China’s CO₂ emission intensity deviates from the long-run equilibrium level; it is corrected by 88.31% in the following year.

For green invention patents (lnGINV), in the long and short terms, green invention patents increased by 1%, China’s CO₂ emission intensity increased by 0.1541%, and decreased by 0.0385% [0.0400 + (−0.0785)] accordingly, and both are significant at the 1% level. It shows that the increase in the number of patents for green inventions, although contributing to the reduction of carbon emission intensity in the short term, inhibits the reduction of carbon emission intensity in the long term, and the long-term effect is stronger than the short-term one. The possible reason is that green invention patents reduce carbon emission intensity in the short term due to their high technology content. However, since the patents for high technology content inventions are usually in the hands of large and high technology enterprises, they may create higher barriers to market access in the long run and, eventually, have a dampening effect on improving carbon emission performance.

On the contrary, green utility model patents (lnGUTY) increased by 1%, leading to a 0.067% increase in China’s CO₂ emission intensity in the short term but a 0.0837% decrease in the long term. The reason for this is that, compared with invention patents, utility model patents are lower in content, and there is a significant patent “bubble” phenomenon, which makes it difficult for them to play a positive role in the short term. Even the negative part is more significant than the positive one. However, due to the strong practicality and easy imitation of green utility model patents, they are conducive to promoting the diffusion and upgrading of green technologies in the long run, thus reducing carbon emission intensity.

The industrial structure (lnIS), measured by the ratio of tertiary to secondary industries, contributes to reducing carbon emission intensity in the short and long terms. However, its long-term effect (−0.3229) is weaker than the short-term effect (−0.7490 + 0.0623 = −0.6967), which indicates that to utilize better the positive effect of industrial restructuring in improving carbon emission performance, it is necessary to further optimize the internal structure of the tertiary industry while promoting the “retreat from two to three.” It indicates that while promoting “shrinking the secondary sector and developing the tertiary sector,” it is necessary to optimize the tertiary sector’s internal structure further to utilize better the positive effect of industrial restructuring in improving carbon emission performance. As the foreign direct investment (lnFDI), measured by FDI dependence, increases (1%), the CO₂ emission intensity increases correspondingly by 0.0529% [0.1088+(−0.0559)] in the short term and by 0.3768% in the long term. It indicates that the larger the share of FDI flows in the GDP of the year, the less favorable it is to improve carbon emission performance. It also reminds us that when introducing foreign investment, the focus should be on setting up more long-term investment attraction programs. Full advantage of the technology spillover is taken from foreign investment introduction rather than just pursuing scale as the goal.

In the short term, as the relative price of energy (lnRP) increases (1%), there is a corresponding decrease in carbon emissions intensity of 0.4641%. In the long run, the strength of its effect is negative (−0.0692) but not significant. It suggests that while raising energy prices is an effective means of reducing carbon emissions, in the long run, there are still some distortions in China’s energy prices, resulting in an energy price mechanism that does not yet play an effective role in allocating resources in the long run. The energy mix measured by the share of non-fossil energy (lnES) fails to play a positive role in effectively reducing carbon emission intensity in the short term (its short-term impact is −0.1026 + 0.1560 = 0.0534), which may be related to China’s coal-dominated resource endowment. Although increasing the share of non-fossil energy helps reducing carbon emissions in the short term, the negative impact on economic development is relatively more significant. It may be related to China’s coal-based resource endowment. However, in the long run, as the share of non-fossil energy increases (1%), the corresponding reduction in carbon intensity −0.4693% plays the strongest role in reducing carbon intensity relative to the other effects examined in this paper. Therefore, increasing the share of non-fossil energy use is crucial for China to achieve the strategic goals of “peak carbon” and “carbon neutrality” in the long run.

When using time series data for cointegration estimation, generally assumed estimated parameters are time-invariant; therefore, further parameter stability tests are required to avoid breakpoints in the time series data and thus affect the stability of the parameters. To this end, this paper applies cumulative sum of recursive residuals (CUSUM) and cumulative sum of squares of recursive residuals (CUSUMQ) based on recursive residuals for parameter stability testing. Figure 1 shows the results. The figure’s upper and lower straight lines are the boundary intervals at a 5% significance level, and the test results indicate that the parameters in the model are stable.

We used FMOLS, DOLS, and CCR to test the robustness of the long-term relationship between variables, and the results are shown in Table 5.

Similar to the long-term results of ARDL, DOLS, FMOLS, and CCR estimates give the same direction of changes in green invention patents, green utility model patents, industrial structure, foreign direct investment, and energy structure as ARDL. Therefore, the results of this study are robust.