Rで産業連関分析

ioanalysis
dibble
R
Japanese
公開

2022年5月31日

産業連関分析は，経済波及効果の算出に広く用いられている分析手法です． 日本では，国や都道府県によって，約5年に1度，産業連関表と呼ばれる統計データが 作成・公開されており，産業連関分析における基礎データとなっています．

これまで，産業連関分析では，Excel・VBAが用いられることが多かったようです．

一方で，近年は，Python・R・Juliaなどのプログラミング言語の普及が進んでいます． これらのプログラミング言語は以下のような特長を持っています．

そのため今後は，産業連関分析においても，これらのプログラミング言語の利用が 進むのではないかと思われます．

ここでは，Rを用いて産業連関分析を行います． Rでは近年，tidyverseなどモダンなデータ分析を行うためのパッケージが 多く提供されており，プログラミング初心者でも習得しやすい言語であると思います．

産業連関表として，e-Statのデータベースで公開されている日本（国）の 2013年・13部門産業連関表を用います． ここで使用するデータは， こちら からダウンロードできます．

産業連関分析の基礎

産業連関分析は一般的に以下の流れに従って行われます．

  1. 産業連関表の整形
  2. 投入係数行列の算出
  3. レオンチェフ逆行列の算出
  4. 経済波及効果の算出

まず，産業連関分析において重要な投入係数行列・レオンチェフ逆行列・経済波及効果の算出方法について解説します．

投入係数行列とは

投入係数は，産業の「クッキングレシピ」として呼ばれており，産業\(j\)の生産物を1単位生産するのに必要な産業\(i\)の生産物の量を表すものです．具体的には，以下のように中間投入\(x_{ij}\)を生産額\(X_j\)（産出額）で割ることで算出できます．

\[ a_{ij}=\frac{x_{ij}}{X_j} \]

産業連関分析では，「クッキングレシピ」に相当する投入係数\(a_{ij}\)に基づく生産額のバランス式（行方向）を連立方程式として解きます．そこで，以下のように，投入係数行列（通常，\(A\)と表される）と呼ばれる行列を作成することで，連立方程式が簡単に解けるようになります．

\[ A = \begin{pmatrix} a_{11} & \cdots & a_{1n} \\ \vdots & \ddots & \vdots \\ a_{n1} & \cdots & a_{nn} \\ \end{pmatrix} \]

レオンチェフ逆行列による経済波及効果の推計について

生産額のバランス式（行方向）は，行列を用いて以下のように表せます．変数の意味は以下の表の通りです．

\[ AX + F + E - M = X \]

変数 意味
\(A\) 投入係数行列
\(X = (X_1, \cdots, X_n) ^ \top\) 生産額ベクトル
\(F = (F_1, \cdots, F_n) ^ \top\) 最終需要ベクトル
\(E = (E_1, \cdots, E_n) ^ \top\) 移輸出ベクトル
\(M = (M_1, \cdots, M_n) ^ \top\) 移輸入ベクトル

経済波及効果の推計では，最終需要の変化が生産額に与える波及効果を算出します．

特に，日本の産業連関表での経済波及効果の推計では，移輸入\(M\)の扱いに注意が必要です（これは，日本表の多くが競争移輸入型表と呼ばれる形式を採用しており，投入係数に移輸入分が含まれているためです）．

最終需要による経済波及効果は，域内の生産額だけでなく域外からの移輸入を誘発すると考えられます．この効果を無視すると経済波及効果を過大評価することにつながるため，通常，投入係数から移輸入相当分を差し引くという処理が行われます．

移輸入は域内需要におおよそ比例すると考えられるため，以下のように，移輸入係数\(\hat{M_i}\)が算出できます．

\[ \hat{M_i} = \frac{M_i}{\sum_{j}a_{ij}X_j + F_i} \]

さらに，行列での計算に適した移輸入係数行列\(\hat{M}\)が，以下のように定義されます．

\[ \hat{M} = \begin{pmatrix} \hat{M_1} & & 0 \\ & \ddots & \\ 0 & & \hat{M_n} \\ \end{pmatrix} \]

以上より，生産額のバランス式（行方向）は，移輸入係数行列\(\hat{M}\)を用いて，以下のように変形されます．ただし，\(I\)は単位行列（対角成分が1，それ以外が0の正方行列）です．

\[ \begin{align} AX + F + E - \hat{M} (AX + F) &= X \\ (I - \hat{M}) (AX + F) + E &= X \end{align} \]

上のバランス式より，経済波及効果の算出式が，以下のように導出されます．ここで，\(\Delta X\)，\(\Delta F\)は，それぞれ，生産額の変化量，最終需要の変化量です．

\[ \begin{align} X &= (I - \hat{M}) (AX + F) + E \\ [I - (I - \hat{M}) A] X &= (I - \hat{M}) F + E \\ X &= [I - (I - \hat{M}) A] ^ {-1} [(I - \hat{M}) F + E] \\ \Delta X &= [I - (I - \hat{M}) A] ^ {-1} (I - \hat{M}) \Delta F \end{align} \]

生産額の変化量\(\Delta X\)の式の右辺の\((I - \hat{M}) \Delta F\)は，最終需要の変化量に自給率\(I - \hat{M}\)を掛けた値となっています．

また，\([I - (I - \hat{M}) A] ^ {-1}\)は，最終需要の変化による直接・間接の波及効果を表す行列であり（開放型または競争移輸入型の）レオンチェフ逆行列と呼ばれています．

以上のように，最終需要の変化量\(\Delta F\)から生産額の変化量\(\Delta X\)を推計するというのが，最も一般的な産業連関分析の方法となっています．

Rによる産業連関分析

産業連関表の整形

ここでは， こちら からダウンロードできる日本の2011年の3部門表（iotable_3sector_2011_wider.csv）を使用します．

こちらの表は，以下のように，日本の2011年の13部門表より作成したもので，単位は「百万円」です．

  • 13部門の産業分類を第1次・第2次・第3次産業に集計（注：「分類不明」を第3次産業に分類）
  • 付加価値部門を1部門に集計
  • 最終需要部門を域内最終需要（finaldemand）・輸出（export）・輸入（import）の3部門に集計

産業連関表のデータ形式は，e-Statのデータベースで提供されている表などを除いて， 行に投入部門（input）・列に産出部門（output）を持つ「横長データ」であることが多いと思われます．

ここでも，以下の通り，まずは横長の産業連関表データを読み込みます．

library(tidyverse)
Warning: package 'ggplot2' was built under R version 4.2.3
Warning: package 'tibble' was built under R version 4.2.3
Warning: package 'dplyr' was built under R version 4.2.3
iotable_wider <- read_csv("iotable_3sector_2011_wider.csv",
                          col_types = cols(.default = "c")) |> 
  
  # input (投入) 列以外を数値に変換
  mutate(across(!input, parse_number))

knitr::kable(iotable_wider)
input industry/01_primary industry/02_secondary industry/03_tertiary finaldemand/04_finaldemand export/05_export import/06_import
industry/01_primary 1456611 7850628 1373767 3869875 47890 -2562809
industry/02_secondary 2715710 161897553 62841827 132924323 54473273 -71673715
industry/03_tertiary 2025270 66811645 155796589 352324555 16423417 -8921553
valueadded/04_valueadded 5838371 106619145 364447740 NA NA NA

データ分析においては，「横長データ」よりも，以下のような「縦長データ」のほうが， 分析しやすい場合が多くあります． ここでも，横長の産業連関表を「縦長データ」に変換します．

iotable <- iotable_wider |>
  
  # input (投入) 列を分類・名称に分割
  separate(input, c("input_type", "input_name"),
           sep = "/") |>
  
  # input (投入) と同様にoutput (産出) の分類・名称列を追加し縦長データに
  pivot_longer(!c(input_type, input_name),
               names_to = c("output_type", "output_name"),
               names_sep = "/",
               values_to = "value_M") |>
  
  # 数値が存在しない行を削除
  drop_na(value_M)

knitr::kable(iotable)
input_type input_name output_type output_name value_M
industry 01_primary industry 01_primary 1456611
industry 01_primary industry 02_secondary 7850628
industry 01_primary industry 03_tertiary 1373767
industry 01_primary finaldemand 04_finaldemand 3869875
industry 01_primary export 05_export 47890
industry 01_primary import 06_import -2562809
industry 02_secondary industry 01_primary 2715710
industry 02_secondary industry 02_secondary 161897553
industry 02_secondary industry 03_tertiary 62841827
industry 02_secondary finaldemand 04_finaldemand 132924323
industry 02_secondary export 05_export 54473273
industry 02_secondary import 06_import -71673715
industry 03_tertiary industry 01_primary 2025270
industry 03_tertiary industry 02_secondary 66811645
industry 03_tertiary industry 03_tertiary 155796589
industry 03_tertiary finaldemand 04_finaldemand 352324555
industry 03_tertiary export 05_export 16423417
industry 03_tertiary import 06_import -8921553
valueadded 04_valueadded industry 01_primary 5838371
valueadded 04_valueadded industry 02_secondary 106619145
valueadded 04_valueadded industry 03_tertiary 364447740

上で構築した表データは，各行のフィルタリングなどが容易にできる一方で， 産業連関分析に用いられる行列計算などに適していません．

そこで，表データの基本的な演算と行列計算を同時に行えるdibbleパッケージを用います． 以下のように，産業連関表をdibbleに変換します．

# pak::pak("UchidaMizuki/dibble")
library(dibble)

iotable <- iotable |>
  dibble_by(input = c(input_type, input_name),
            output = c(output_type, output_name),
            
            # "_"で列名を分割してinput (投入)・output (産出) 軸を設定
            .names_sep = "_")

iotable
# A dibble:   24 x 1
# Dimensions: input [4], output [6]
# Measures:   value_M
   input$type $name        output$type $name            value_M
   <chr>      <chr>        <chr>       <chr>              <dbl>
 1 industry   01_primary   industry    01_primary       1456611
 2 industry   01_primary   industry    02_secondary     7850628
 3 industry   01_primary   industry    03_tertiary      1373767
 4 industry   01_primary   finaldemand 04_finaldemand   3869875
 5 industry   01_primary   export      05_export          47890
 6 industry   01_primary   import      06_import       -2562809
 7 industry   02_secondary industry    01_primary       2715710
 8 industry   02_secondary industry    02_secondary   161897553
 9 industry   02_secondary industry    03_tertiary     62841827
10 industry   02_secondary finaldemand 04_finaldemand 132924323
# ℹ 14 more rows

投入係数行列の算出

産業の「クッキングレシピ」と呼ばれる投入係数行列\(A\)を以下のように中間投入を生産額で割って算出します．

注：dibbleではブロードキャストが自動で行われますが，安全のため，ブロードキャストを行う際に，警告を発するように設計されています．そのため，broadcast()でブロードキャスト後の軸名c("input", "output")を与えて警告が出ないようにする必要があります．

# 生産額
total_input <- iotable |>
  filter(output$type == "industry") |>
  apply("output", sum)

# 中間投入
interindustry <- iotable |>
  filter(input$type == "industry",
         output$type == "industry")

# 投入係数
inputcoeff <- broadcast(interindustry / total_input,
                        c("input", "output"))

inputcoeff
# A dibble:   9
# Dimensions: input [3], output [3]
  input$type $name        output$type $name              .
  <chr>      <chr>        <chr>       <chr>          <dbl>
1 industry   01_primary   industry    01_primary   0.121  
2 industry   01_primary   industry    02_secondary 0.0229 
3 industry   01_primary   industry    03_tertiary  0.00235
4 industry   02_secondary industry    01_primary   0.226  
5 industry   02_secondary industry    02_secondary 0.472  
6 industry   02_secondary industry    03_tertiary  0.108  
7 industry   03_tertiary  industry    01_primary   0.168  
8 industry   03_tertiary  industry    02_secondary 0.195  
9 industry   03_tertiary  industry    03_tertiary  0.267

レオンチェフ逆行列の算出

経済波及効果を表すレオンチェフ逆行列は以下のように，移輸入係数と投入係数を用いて算出できます．

注：solve()で逆行列を算出すると行列の軸名が入れ替わるため注意してください．

# 域内需要
localdemand <- iotable |>
  filter(input$type == "industry",
         !output$type %in% c("export", "import")) |>
  apply("input", sum)

# (移) 輸入
import <- iotable |>
  filter(input$type == "industry",
         output$type == "import") |>
  apply("input", sum)
# 符号を正に
import <- -import

# (移) 輸入係数
importcoeff <- import / localdemand

I <- eye(inputcoeff) # 単位行列
M <- importcoeff     # 移輸入係数ベクトル (broadcastが行われるため行列でなくてよい)
A <- inputcoeff      # 投入係数行列

# レオンチェフ逆行列
leontiefinv <- broadcast(I - (1 - M) * A,
                         c("input", "output")) |>
  solve()

leontiefinv
# A dibble:   9
# Dimensions: output [3], input [3]
  output$type $name        input$type $name              .
  <chr>       <chr>        <chr>      <chr>          <dbl>
1 industry    01_primary   industry   01_primary   1.12   
2 industry    01_primary   industry   02_secondary 0.0361 
3 industry    01_primary   industry   03_tertiary  0.00716
4 industry    02_secondary industry   01_primary   0.374  
5 industry    02_secondary industry   02_secondary 1.68   
6 industry    02_secondary industry   03_tertiary  0.197  
7 industry    03_tertiary  industry   01_primary   0.348  
8 industry    03_tertiary  industry   02_secondary 0.445  
9 industry    03_tertiary  industry   03_tertiary  1.41

経済波及効果の算出

こちら からダウンロードできる最終需要がそれぞれ百万円ずつ増加する（finaldemand_change_3sector.csv）ケースで経済波及効果を算出しています．

# 最終需要変化量
finaldemand_change <- read_csv("finaldemand_change_3sector.csv",
                               col_types = cols(.default = "c",
                                                value_M = "n")) |> 
  dibble_by(input = c(input_type, input_name),
            .names_sep = "_")

L <- leontiefinv         # レオンチェフ逆行列
M <- importcoeff         # 移輸入係数
FD <- finaldemand_change # 最終需要変化量

# 経済波及効果
spillover <- L %*% ((1 - M) * FD)

spillover
# A dibble:   3
# Dimensions: output [3]
  output$type $name            .
  <chr>       <chr>        <dbl>
1 industry    01_primary   0.958
2 industry    02_secondary 1.85 
3 industry    03_tertiary  2.03 
theme_set(theme_light())

spillover |> 
  as_tibble(n = "value_M") |> 
  unpack(output, 
         names_sep = "_") |> 
  ggplot(aes(output_name, value_M,
             fill = output_name)) +
  geom_col(show.legend = FALSE) +
  scale_fill_brewer(palette = "Set2")

まとめ

Rを用いた産業連関分析の方法について紹介しました．

ここまでの計算を，jpioにパッケージ形式でまとめました．以下のように，ここまでの計算と同様の計算を行うことができます．

# pak::pak("UchidaMizuki/jpio")

# 産業連関表
iotable <- read_csv("iotable_3sector_2011.csv",
                    col_types = cols(.default = "c",
                                     value_M = "n")) |> 
  jpio::as_iotable()

iotable
# A dibble:   24
# Dimensions: input [4], output [6]
   input$type $name        output$type $name                  .
   <fct>      <chr>        <fct>       <chr>              <dbl>
 1 industry   01_primary   industry    01_primary       1456611
 2 industry   01_primary   industry    02_secondary     7850628
 3 industry   01_primary   industry    03_tertiary      1373767
 4 industry   01_primary   finaldemand 04_finaldemand   3869875
 5 industry   01_primary   export      05_export          47890
 6 industry   01_primary   import      06_import       -2562809
 7 industry   02_secondary industry    01_primary       2715710
 8 industry   02_secondary industry    02_secondary   161897553
 9 industry   02_secondary industry    03_tertiary     62841827
10 industry   02_secondary finaldemand 04_finaldemand 132924323
# ℹ 14 more rows
# 投入係数
jpio::input_coef(iotable)
# A dibble:   9
# Dimensions: input [3], output [3]
  input$type $name        output$type $name              .
  <fct>      <chr>        <fct>       <chr>          <dbl>
1 industry   01_primary   industry    01_primary   0.121  
2 industry   01_primary   industry    02_secondary 0.0229 
3 industry   01_primary   industry    03_tertiary  0.00235
4 industry   02_secondary industry    01_primary   0.226  
5 industry   02_secondary industry    02_secondary 0.472  
6 industry   02_secondary industry    03_tertiary  0.108  
7 industry   03_tertiary  industry    01_primary   0.168  
8 industry   03_tertiary  industry    02_secondary 0.195  
9 industry   03_tertiary  industry    03_tertiary  0.267  
# レオンチェフ逆行列
jpio::leontief_inv(iotable)
# A dibble:   9
# Dimensions: output [3], input [3]
  output$type $name        input$type $name              .
  <fct>       <chr>        <fct>      <chr>          <dbl>
1 industry    01_primary   industry   01_primary   1.12   
2 industry    01_primary   industry   02_secondary 0.0361 
3 industry    01_primary   industry   03_tertiary  0.00716
4 industry    02_secondary industry   01_primary   0.374  
5 industry    02_secondary industry   02_secondary 1.68   
6 industry    02_secondary industry   03_tertiary  0.197  
7 industry    03_tertiary  industry   01_primary   0.348  
8 industry    03_tertiary  industry   02_secondary 0.445  
9 industry    03_tertiary  industry   03_tertiary  1.41   
# 経済波及効果
jpio::spillover_effect(iotable,
                       list(`01_primary` = 1,
                            `02_secondary` = 1,
                            `03_tertiary` = 1))
# A dibble:   3
# Dimensions: output [3]
  output$type $name            .
  <fct>       <chr>        <dbl>
1 industry    01_primary   0.958
2 industry    02_secondary 1.85 
3 industry    03_tertiary  2.03 
# スカイラインチャート
jpio::skyline_chart(iotable, 
                    ylim = c(-0.5, NA))