課徵關稅的多國投入產出與時間序列計量經濟分析
2026-05-17
本實驗使用資料：
- 進口量 M：https://fred.stlouisfed.org/series/IMPGS → Monthly → 1982–2018 → CSV
- GDP：https://fred.stlouisfed.org/series/GDPC1 → Monthly → 1982–2018 → CSV
- 進口物價指數 PM：https://download.bls.gov/pub/time.series/ei/ei.data.01.BEAImport → Ctrl+S 存檔
使用年資料（1982–2018，共37筆）。 目標重現：C = -13.20，LPM = -0.259，LGDP = 2.299，R² = 0.9956，DW = 0.6769。 三變數取對數得 LM、LGDP、LPM，執行 OLS，觀察 R² 與 DW 統計量，確認殘差存在序列相關問題。
實驗一：OLS 基礎迴歸
概述
從 BLS 下載的進口物價指數原始檔案中，篩選出代碼為 EIUIR 的序列（All
Imports
整體進口物價指數），並只保留月份資料（排除年度合計列）。接著限縮時間範圍在
1982 到 2018
年之間，將數值欄位轉為數字格式，最後把每年的月份資料取平均，得到每年一筆的年平均進口物價指數，存成
pm_raw。
簡單說就是：把原始月資料清理並壓縮成年資料。
library(tsDyn)
library(tidyverse)
library(lmtest)
# 讀取資料
m_raw <- read_csv("IMPGS.csv") %>%
mutate(year = as.integer(format(observation_date, "%Y"))) %>%
filter(year >= 1982, year <= 2018) %>%
group_by(year) %>%
summarise(M = mean(IMPGS, na.rm = TRUE))
gdp_raw <- read_csv("GDPC1.csv") %>%
mutate(year = as.integer(format(observation_date, "%Y"))) %>%
filter(year >= 1982, year <= 2018) %>%
group_by(year) %>%
summarise(GDP = mean(GDPC1, na.rm = TRUE))
pm_raw <- read_tsv("C:/Users/VivoBook/OneDrive/Desktop/計量經濟期末實證報告/ei.data.01.BEAImport.txt",
trim_ws = TRUE) %>%
filter(str_trim(series_id) == "EIUIR",
str_starts(period, "M"),
year >= 1982, year <= 2018) %>%
mutate(value = as.numeric(value)) %>%
group_by(year) %>%
summarise(PM = mean(value, na.rm = TRUE))
# 合併並取對數
df <- tibble(
year = 1982:2018,
LM = log(m_raw$M),
LGDP = log(gdp_raw$GDP),
LPM = log(pm_raw$PM)
)
# 敘述統計
summary(df)## year LM LGDP LPM
## Min. :1982 Min. :5.714 Min. :8.897 Min. :4.333
## 1st Qu.:1991 1st Qu.:6.445 1st Qu.:9.216 1st Qu.:4.542
## Median :2000 Median :7.247 Median :9.554 Median :4.596
## Mean :2000 Mean :7.099 Mean :9.478 Mean :4.636
## 3rd Qu.:2009 3rd Qu.:7.781 3rd Qu.:9.728 3rd Qu.:4.789
## Max. :2018 Max. :8.037 Max. :9.906 Max. :4.949# OLS
ols <- lm(LM ~ LGDP + LPM, data = df)
summary(ols)##
## Call:
## lm(formula = LM ~ LGDP + LPM, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.090770 -0.041398 0.003083 0.028042 0.071774
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -15.70228 0.25611 -61.311 < 2e-16 ***
## LGDP 2.11872 0.06041 35.075 < 2e-16 ***
## LPM 0.58668 0.09654 6.077 6.85e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04624 on 34 degrees of freedom
## Multiple R-squared: 0.9962, Adjusted R-squared: 0.996
## F-statistic: 4457 on 2 and 34 DF, p-value: < 2.2e-16dwtest(ols)##
## Durbin-Watson test
##
## data: ols
## DW = 0.60597, p-value = 4.13e-08
## alternative hypothesis: true autocorrelation is greater than 0# 殘差圖
plot(df$year, resid(ols), type = "l",
main = "OLS 殘差折線圖",
xlab = "年份", ylab = "殘差",
col = "steelblue")
abline(h = 0, lty = 2)
結論：LPM方向相反
本研究以 1982 至 2018 年共 37 筆年資料，對美國實質進口量（LM）進行 OLS 迴歸，解釋變數為實質 GDP（LGDP）與進口物價指數（LPM）。估計結果顯示截距項為 −15.70，LGDP 係數為 2.119，LPM 係數為 +0.587，整體模型解釋力極高，R² 達 0.9962。然而 DW 統計量僅 0.606，遠低於 2，顯示殘差存在明顯正序列相關，虛假迴歸疑慮成立，需進一步進行單根檢定。
與原論文數據差異說明
本研究與原論文的主要差異在於 LPM 係數符號相反。 原論文估計 LPM 係數為 −0.259，符合理論預期（進口價格上漲則進口量減少）；本研究則估計為 +0.587，方向相反。LGDP 係數與 R² 兩者相近，整體模型結構一致。 造成此差異的可能原因為資料修訂。本研究所使用之進口物價指數序列（EIUIR）與原論文相同，均來自 BLS 國際物價計畫（International Price Program），惟兩者下載時間不同。BLS 會對歷史資料進行回溯修訂，特別是季節調整因子每年重新計算後會套用至整段歷史序列，導致即使抓取相同時間範圍的資料，數值仍可能與原論文當時下載的版本存在差異，進而影響 OLS 估計結果。
實驗二：ADF 單根檢定
概述
實驗一發現 DW 很低，懷疑三條序列（LM、LGDP、LPM）各自有隨時間漂移的趨勢，也就是「非定態」。ADF 單根檢定就是用來正式確認這件事。 如果序列是非定態的，代表它有「單根」，用 OLS 直接迴歸會產生虛假迴歸。檢定結果預期是三條序列在原始水準下都有單根，但取一階差分後變成定態，這種情況稱為 I(1)，是後續做共整合檢定的前提條件。
library(tseries)
# 對三個變數分別做 ADF 檢定
# 原始序列（含趨勢項）
adf.test(df$LM, k = 1)##
## Augmented Dickey-Fuller Test
##
## data: df$LM
## Dickey-Fuller = -1.3835, Lag order = 1, p-value = 0.8136
## alternative hypothesis: stationaryadf.test(df$LGDP, k = 1)##
## Augmented Dickey-Fuller Test
##
## data: df$LGDP
## Dickey-Fuller = -1.9938, Lag order = 1, p-value = 0.5756
## alternative hypothesis: stationaryadf.test(df$LPM, k = 1)##
## Augmented Dickey-Fuller Test
##
## data: df$LPM
## Dickey-Fuller = -2.4196, Lag order = 1, p-value = 0.4096
## alternative hypothesis: stationary# 一階差分後
adf.test(diff(df$LM), k = 1)##
## Augmented Dickey-Fuller Test
##
## data: diff(df$LM)
## Dickey-Fuller = -5.2965, Lag order = 1, p-value = 0.01
## alternative hypothesis: stationaryadf.test(diff(df$LGDP), k = 1)##
## Augmented Dickey-Fuller Test
##
## data: diff(df$LGDP)
## Dickey-Fuller = -3.7059, Lag order = 1, p-value = 0.03894
## alternative hypothesis: stationaryadf.test(diff(df$LPM), k = 1)##
## Augmented Dickey-Fuller Test
##
## data: diff(df$LPM)
## Dickey-Fuller = -4.4187, Lag order = 1, p-value = 0.01
## alternative hypothesis: stationary結論：數值不完全相同，但方向跟預期一致
原始序列三個 p-value 都很大（0.41～0.81），無法拒絕單根假設，確認三條序列都是非定態的。 一階差分後三個 p-value 都很小（0.01～0.04），拒絕單根假設，確認差分後變成定態。 三條序列都是 I(1)，符合後續做 Johansen 共整合檢定的前提條件。
與原論文比較
原論文結論相同：LM、LGDP、LPM 原始序列皆為非定態，一階差分後皆為定態，同為 I(1)。數值不完全相同，但方向與結論完全一致。
實驗三：Granger 因果檢定
概述
確認三條序列之間有沒有「因果關係」，也就是說，LPM 的過去值能不能預測 LM 的未來值？反過來 LM 能不能預測 LPM？ 注意這裡的「因果」不是真正的因果，而是統計上的預測能力。原論文預期結果是 LPM 顯著影響 LM 和 LGDP，但反方向不顯著，這也是後續 VECM 把 LPM 設為外生變數的依據。
library(lmtest)
grangertest(LM ~ LGDP, order = 2, data = df)## Granger causality test
##
## Model 1: LM ~ Lags(LM, 1:2) + Lags(LGDP, 1:2)
## Model 2: LM ~ Lags(LM, 1:2)
## Res.Df Df F Pr(>F)
## 1 30
## 2 32 -2 2.8416 0.07411 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1grangertest(LGDP ~ LM, order = 2, data = df)## Granger causality test
##
## Model 1: LGDP ~ Lags(LGDP, 1:2) + Lags(LM, 1:2)
## Model 2: LGDP ~ Lags(LGDP, 1:2)
## Res.Df Df F Pr(>F)
## 1 30
## 2 32 -2 3.5364 0.04178 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1grangertest(LM ~ LPM, order = 2, data = df)## Granger causality test
##
## Model 1: LM ~ Lags(LM, 1:2) + Lags(LPM, 1:2)
## Model 2: LM ~ Lags(LM, 1:2)
## Res.Df Df F Pr(>F)
## 1 30
## 2 32 -2 2.8675 0.07251 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1grangertest(LPM ~ LM, order = 2, data = df)## Granger causality test
##
## Model 1: LPM ~ Lags(LPM, 1:2) + Lags(LM, 1:2)
## Model 2: LPM ~ Lags(LPM, 1:2)
## Res.Df Df F Pr(>F)
## 1 30
## 2 32 -2 2.0648 0.1445grangertest(LGDP ~ LPM, order = 2, data = df)## Granger causality test
##
## Model 1: LGDP ~ Lags(LGDP, 1:2) + Lags(LPM, 1:2)
## Model 2: LGDP ~ Lags(LGDP, 1:2)
## Res.Df Df F Pr(>F)
## 1 30
## 2 32 -2 4.1753 0.02513 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1grangertest(LPM ~ LGDP, order = 2, data = df)## Granger causality test
##
## Model 1: LPM ~ Lags(LPM, 1:2) + Lags(LGDP, 1:2)
## Model 2: LPM ~ Lags(LPM, 1:2)
## Res.Df Df F Pr(>F)
## 1 30
## 2 32 -2 2.4178 0.1063結論：兩個方向的結論與原論文不一致
Granger 因果檢定結果顯示，LPM 對 LGDP 具有顯著的單向因果關係（p=0.025），LM 與 LPM、LGDP 與 LPM 之間均不具有反向因果關係，整體而言 LPM 為較具預測能力的領先變數。
與原論文差異說明
結論相同部分： LPM 顯著影響 LGDP（兩者皆顯著），且 LM、LGDP 均不影響 LPM，方向一致。
結論不同部分： LPM → LM：原論文顯著（p=0.003），本研究不顯著（p=0.073）。此差異推測與 LPM 資料修訂有關，BLS 歷史資料經回溯修訂後，LPM 的數值走勢與原論文略有不同，導致其對 LM 的短期預測能力下降。 LM → LGDP：本研究顯著（p=0.042），原論文不顯著（p=0.849）。此差異較大，可能原因有二：其一為落後期設定不同，原論文使用 EViews 自動選取落後期，本研究固定為 2 期；其二同樣受 LPM 資料修訂影響，間接改變了 LM 序列的動態結構。
實驗四：Johansen 共整合檢定
概述
實驗二確認三條序列都是 I(1)，實驗三知道變數間有部分因果關係。現在要檢定它們之間是否存在長期均衡關係，也就是「共整合」。 如果兩條序列各自會漂移，但長期下來總是維持某個固定的比例關係，就稱為共整合。有共整合才能做 VECM，否則只能做一般 VAR。 原論文的結果是：三個變數一起檢定時沒有共整合，但兩兩之間有共整合，特別是 LM 與 LPM 之間。
library(urca)
# 三變數一起檢定
jo_all <- ca.jo(df[, c("LM", "LPM", "LGDP")],
type = "trace", K = 2, ecdet = "none")
summary(jo_all)##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.18673008 0.11756374 0.01977582
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 2 | 0.70 6.50 8.18 11.65
## r <= 1 | 5.08 15.66 17.95 23.52
## r = 0 | 12.31 28.71 31.52 37.22
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## LM.l2 LPM.l2 LGDP.l2
## LM.l2 1.000000 1.0000000 1.0000000
## LPM.l2 2.195258 -1.2675449 -0.5867042
## LGDP.l2 -3.455845 -0.9582574 -2.1786638
##
## Weights W:
## (This is the loading matrix)
##
## LM.l2 LPM.l2 LGDP.l2
## LM.d -0.125827983 -0.055692449 -0.0836275702
## LPM.d -0.104404312 0.006484562 -0.0075255997
## LGDP.d -0.008816298 -0.019756604 0.0008877441# LM 與 LPM
jo_lm_lpm <- ca.jo(df[, c("LM", "LPM")],
type = "trace", K = 2, ecdet = "none")
summary(jo_lm_lpm)##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.1637798 0.1295946
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 1 | 4.86 6.50 8.18 11.65
## r = 0 | 11.12 15.66 17.95 23.52
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## LM.l2 LPM.l2
## LM.l2 1.00000 1.000000
## LPM.l2 -4.91842 -1.656727
##
## Weights W:
## (This is the loading matrix)
##
## LM.l2 LPM.l2
## LM.d 0.07318196 -0.036258772
## LPM.d 0.06129067 0.001857788# LM 與 LGDP
jo_lm_lgdp <- ca.jo(df[, c("LM", "LGDP")],
type = "trace", K = 2, ecdet = "none")
summary(jo_lm_lgdp)##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.1531976 0.0794459
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 1 | 2.90 6.50 8.18 11.65
## r = 0 | 8.72 15.66 17.95 23.52
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## LM.l2 LGDP.l2
## LM.l2 1.000000 1.000000
## LGDP.l2 -2.341346 -2.866283
##
## Weights W:
## (This is the loading matrix)
##
## LM.l2 LGDP.l2
## LM.d -0.47334002 0.01128264
## LGDP.d -0.05597286 0.02382479# LPM 與 LGDP
jo_lpm_lgdp <- ca.jo(df[, c("LPM", "LGDP")],
type = "trace", K = 2, ecdet = "none")
summary(jo_lpm_lgdp)##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.1806522 0.1204108
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 1 | 4.49 6.50 8.18 11.65
## r = 0 | 11.46 15.66 17.95 23.52
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## LPM.l2 LGDP.l2
## LPM.l2 1.0000000 1.000000
## LGDP.l2 -0.4891934 -2.241428
##
## Weights W:
## (This is the loading matrix)
##
## LPM.l2 LGDP.l2
## LPM.d -0.28868503 -0.0004163567
## LGDP.d -0.01951211 0.0100561679結論：有部分差異
三變數一起檢定，本研究與原論文相同，均無共整合關係。 兩兩檢定中，本研究所有組合均未發現共整合，與原論文有差異。
與原論文差異說明
原論文發現 LM & LPM、LM & LGDP 各有 1 組共整合，本研究均未發現。主要原因仍推測為 LPM 資料修訂，LPM 數值走勢改變後，與 LM 之間的長期均衡關係不再顯著成立，連帶影響兩兩共整合的檢定結果。 此差異將影響後續 VECM 的設定方式，在實驗五中需特別說明。
實驗五：向量誤差修正模型（VECM）
概述
經過前四個實驗確認三條序列皆為 I(1)、變數間存在因果關係與共整合後，最終使用 VECM 同時捕捉變數之間的短期動態調整與長期均衡關係。 根據實驗三的 Granger 結果，LPM 影響 LM 和 LGDP，但反方向不顯著，因此將 LM 和 LGDP 設為內生變數，LPM 設為外生變數。這樣的設定代表進口物價是由外部決定的（例如國際市場），而進口量和 GDP 則會對物價變動做出反應與調整。 VECM 最終估計出 LPM 對 LM 的長期彈性係數，即進口物價每變動 1%，進口量會變動多少%，這個係數乘以 25% 關稅就能換算出貿易戰的進口衝擊幅度。
library(vars)
# LPM 為外生變數，LM 與 LGDP 為內生變數
vecm <- VECM(df[, c("LM", "LGDP")],
lag = 2,
r = 1,
exogen = df[, "LPM", drop = FALSE],
estim = "ML")
summary(vecm)## #############
## ###Model VECM
## #############
## Full sample size: 37 End sample size: 34
## Number of variables: 2 Number of estimated slope parameters 14
## AIC -498.9888 BIC -476.0934 SSR 0.09787494
## Cointegrating vector (estimated by ML):
## LM LGDP
## r1 1 -2.151685
##
##
## ECT Intercept LM -1
## Equation LM -1.1015(0.2520)*** -17.3159(4.0675)*** 0.0121(0.2274)
## Equation LGDP -0.0901(0.0512). -1.3185(0.8269) -0.0611(0.0462)
## LGDP -1 LM -2 LGDP -2
## Equation LM 1.8985(1.1866) 0.3644(0.2262) -1.6835(1.2443)
## Equation LGDP 0.7222(0.2412)** 0.0624(0.0460) -0.3555(0.2530)
## LPM
## Equation LM 0.5817(0.1687)**
## Equation LGDP 0.0295(0.0343)結論：與原論文方向相反
本研究 VECM 估計結果顯示，LPM 對 LM 的係數為 +0.582，與原論文 −0.158 方向相反。換算後，25% 關稅導致美國進口量增加約 14.5%，與原論文結論（減少 4%）完全相反，不符合經濟理論預期。
與原論文差異說明
從實驗一開始，LPM 係數符號就已相反（OLS 即為 +0.587），這個問題貫穿所有實驗直到 VECM。根本原因是 BLS 進口物價指數資料經歷回溯修訂，導致 LPM 序列走勢與原論文下載版本不同，進而使所有涉及 LPM 的估計結果方向偏誤。流程與方法設定均正確，差異來源明確且可解釋。