2015年に書かれたRを解読する

背景

筆者は学生の時から同じPCを使っているため、PCを漁ると太古の昔の文法で書かれたRコードを発掘できる。
その中に、先輩に書いて頂いたものかつ、内容が全くわからないコードを見つけたため、解読ついでに現在の文法へのupdateを行うことにした。

課題

• 大昔のコードを解読し現代の文法にupdateすることで、Rへの理解を深める
• Rのここ10年での進化を実感し感謝する

本題

元の2015年版

RAS <- function( A, x, z, y, eps=1.0e-8, itMax=500, verbose=F ){
  itCount  <- 0
  converge <- F
  invalid <- F

  length_X <- length(x)
  ncol_A <- ncol(A)
  if(length_X != ncol_A) {invalid <- T}

   if ( length( z ) != ncol( A ) ) invalid <- T
  if ( length( y ) != nrow( A ) ) invalid <- T

  
 
  if ( invalid ) {
    cat( "Data is not consistent:¥n" )
    cat( "  A: ", nrow( A ), "x", ncol( A ), "¥n" )
    cat( "  x: 1 x", length( x ), "¥n" )
    cat( "  z: 1 x", length( z ), "¥n" )
    cat( "  y: ", length( y ), "x 1¥n" )
    return( "converge"=converge, "iteration"=itCount, "Table"=0 )
  }

  X <- diag( x )

  while ( !converge ) {
    itCount <- itCount + 1
    if ( itCount > itMax ) break
    B <- A %*% X
    s <- z / apply( B, 2, sum )
    s[ is.na( s ) | is.nan( s ) | is.infinite( s ) ] <- 0
    A <- A %*% diag( s )
    B <- A %*% X
    r <- y / apply( B, 1, sum )
    r[ is.na( r ) | is.nan( r ) | is.infinite( r ) ] <- 0
    A <- diag( r ) %*% A
    rMax <- max( abs( r ) ) - 1
    sMax <- max( abs( s ) ) - 1
    if ( verbose ) {
      cat( "iteration: ", itCount, "¥n" )
      cat( "s¥n" )
      print( s )
      cat( "r¥n" )
      print( r )
    }
    if ( rMax <= eps & sMax <= eps ) converge <- T
  }
  cat( "Total Iteration: ", itCount, "¥n" )
  cat( "Convergence: ", converge, "¥n" )
  cat( "s¥n" )
  print( s )
  cat( "r¥n" )
  print( r )
  return( "converge"=converge, "iteration"=itCount, "Table"=A %*% X )
}

RAS(a,x,z,y)
RAS(A=a,x=x,z=z,y=y)

RAS <- function( A, x, z, y, eps=1.0e-8, itMax=500, verbose=FALSE ){
  itCount  <- 0
  converge <- FALSE
  invalid <- FALSE

length_X <- length(x)
ncol_A <- ncol(A)

  if ( length( x ) != ncol( A ) ) invalid <- T
  if ( length( z ) != ncol( A ) ) invalid <- T
  if ( length( y ) != nrow( A ) ) invalid <- T

  return(invalid)
}

このRのコードは、 「RAS法」を実装している関数であった。RAS法は、与えられた行列の行と列の合計が特定の目標ベクトル（xとy）に一致するように行列を調整するために使用される。この手法は、計量経済学界隈で産業連関表(詳しくは後述)の作成に使われる。要は、そこそこマニアックな手法である。

コードの概要は以下の通り：

1. 初期化：反復回数カウンター(itCount)、収束フラグ(converge)、無効フラグ(invalid)を初期化。
2. 入力チェック：入力された行列Aとベクトルx、z、yの長さが一致するかどうかをチェックし、一致しない場合はエラーメッセージを出力して関数を終了。
3. RASアルゴリズムの実行：行列Aに対してRASアルゴリズムを適用し、行と列の合計がベクトルxとyに一致するように調整。この過程は、収束するか最大反復回数に達するまで続けられる。
4. 収束チェック：各反復で、行と列のスケーリングファクターが十分に小さいかどうかをチェックし、小さい場合は収束したとみなす。
5. 結果の出力：収束したかどうか、反復回数、および調整された行列を出力する。

2024年版

tidyverseやdplyrは主にデータフレーム操作に使用されるため、行列操作にはあまり適していない。したがって、ここでは基本的なR関数を使用してコードをupdateした。

# RAS法による行列バランス調整
# A: バランスを取る必要がある入力行列
# x: 目標とする列の合計
# z: 目標とする行の合計
# y: 目標とする行の合計（zと重複している可能性があるため、エラーの可能性あり）
# eps: 収束基準
# itMax: 最大反復回数
# verbose: 反復の詳細を出力するフラグ
RAS <- function(A, x, z, y, eps = 1.0e-8, itMax = 500, verbose = FALSE) {
  itCount <- 0  # 反復カウンタの初期化
  converge <- FALSE  # 収束フラグの初期化
  invalid <- FALSE   # 有効性フラグの初期化

  # 入力の次元が一致しているかチェック
  if (length(x) != ncol(A) || length(z) != ncol(A) || length(y) != nrow(A)) {
    invalid <- TRUE
  }

  # データが一致しない場合はエラーメッセージを出力し、終了
  if (invalid) {
    cat("データが一致しません:\n")
    cat("  A: ", nrow(A), "x", ncol(A), "\n")
    cat("  x: 1 x", length(x), "\n")
    cat("  z: 1 x", length(z), "\n")
    cat("  y: ", length(y), "x 1\n")
    return(list(converge = converge, iteration = itCount, Table = NULL))
  }

  # 目標列合計で対角行列Xを初期化
  X <- diag(x)

  # RAS手順の反復
  while (!converge && itCount <= itMax) {
    itCount <- itCount + 1
    B <- A %*% X  # 行列Bの更新
    s <- z / colSums(B)  # 列のスケーリング係数を計算
    s[is.na(s) | is.nan(s) | is.infinite(s)] <- 0  # sのNA, NaN, Infを0に置換
    A <- A %*% diag(s)  # Aの列をスケール
    B <- A %*% X  # 行列Bの更新
    r <- y / rowSums(B)  # 行のスケーリング係数を計算
    r[is.na(r) | is.nan(r) | is.infinite(r)] <- 0  # rのNA, NaN, Infを0に置換
    A <- diag(r) %*% A  # Aの行をスケール
    # 収束チェック
    rMax <- max(abs(r) - 1)
    sMax <- max(abs(s) - 1)
    if (verbose) {
      cat("反復: ", itCount, "\n")
      cat("s\n")
      print(s)
      cat("r\n")
      print(r)
    }
    if (rMax <= eps && sMax <= eps) converge <- TRUE
  }
  # 最終的な詳細を出力
  cat("総反復回数: ", itCount, "\n")
  cat("収束: ", converge, "\n")
  if (verbose) {
    cat("s\n")
    print(s)
    cat("r\n")
    print(r)
  }
  # 結果をリストとして返す
  return(list(converge = converge, iteration = itCount, Table = A %*% X))
}

修正内容

• apply 関数の代わりに colSums と rowSums を使用してみた。
• return ステートメントで NULL を返すようにした。
• & 演算子の代わりに && を使用した。
• verbose フラグが TRUE の場合にのみ、最終的な s と r を出力するようにした。

2015年版が生み出された経緯

筆者は当時学部3年生で、計量経済学分野のゼミで産業連関表をあーだこーだしていた。
係数や逆行列の計算はExcelでできたものの、反復のプロセスを要するRAS法の部分についてはExcelでは処理できなかったため、当時既につよつよプログラミングマンだった先輩にコードを書いて頂いた。

産業連関表・産業連関分析とは

産業連関表・産業連関分析は、計量経済学界隈の中でもハードコアというか、2024年現在めちゃくちゃホットな分野ではないのではないかと推測するため、簡単に概要を記載する。

産業連関表は、ある地域内の一定期間（通常は一年間）の経済活動・産業間取引を一つのマトリックスに示した統計表であり、地域内で生産される商品の投入構造と需要構造とを両面から捉えることで、生産部門間の相互依存の姿を体系的に記述している。産業連関表によって、地域内の生産額や付加価値額や最終需要額などの作成対象年次の地域経済構造を把握できるだけでなく、産業連関表の構造を利用して経済構造の将来予測、価格分析、政策・事業の経済効果の測定などの産業連関分析をすることも可能である。

産業連関分析を行ううえで基本になるのは取引基本表・投入係数表・逆行列係数表の 3 つである。

■取引基本表

各産業相互間や家計・政府などの消費者との間で財・サービスがどのように生産され、販売されたかを示したものである。この表を横（行）に読めば販路構成が分かり、縦（列）に読めば投入した費用の構成を知ることができる。なお、取引基本表は各部門全て、縦に見た投入額合計と横に見た産出額合計が一致するように作成されている。

■投入係数表

ある産業で一単位の生産を行う際に必要となる原材料等の投入量の割合である投入係数を一覧表にして示したものである。

■逆行列係数表

ある産業に対する最終需要が一単位発生した時、その需要を満たすために必要な産業別の生産が最終的にどれくらいになるかという、生産波及の大きさを示す係数を一覧表にして示したものである。

産業連関分析 推しコメント：

最近のトピックで例を上げるなら、「令和6年能登半島地震の経済波及効果」「2025大阪・関西万博の経済波及効果」といったたぐいの分析が可能である。
「経済波及効果」として直接の一次効果のみならずn次効果まで算出可能であり、地域間の分析も可能である。更には、CO2排出量等の環境分析にも応用可能である。
産業連関表は、「1枚の平面に全世界の経済(スケールは地球でも国でも都道府県でもOK)を表現できているだけでおもろいのに、しかも色々分析して遊べる！」という、ダイナミックで愉快なツールです！

おまけ：プログラミングにおける「function」とは

2015年版の作者に確認したところ、当時筆者が「そもそもfunctionってなんすか？」とぬかしていたとのことだったので(アホの子すぎる)、改めて整理しておく。

プログラミングにおける「function」（関数）は、特定のタスクを実行するためのコードのブロックである。関数は、一連の命令をカプセル化し、それを繰り返し使用することを可能にします。関数は、入力（引数やパラメータと呼ばれる）を受け取り、処理を行い、結果（戻り値）を返すことができる。

関数の主な利点は次のとおり：

1. 再利用性: 同じコードを何度も書く代わりに、関数を定義して必要なときに呼び出すことができる。
2. モジュール性: 大きなプログラムを小さな、管理しやすい部分に分割することができる。
3. メンテナンス: 関数内のコードを変更すると、その関数を使用するプログラム全体に変更が適用される。

関数の例：
① python

def shinjiro():
    if True:
        return True
print(shinjiro())

② R

shinjiro <- function() {
  if (TRUE) {
    return(TRUE)
  }
}
print(shinjiro())

所感

今回のコードはデータフレーム操作をするものではないため、この10年で生まれたtidyverseやdplyrの素晴らしさを証明できる内容にはならなかった。

2024年に産業連関分析をしている人はもう絶滅危惧種かもしれませんが、今どこかで産業連関分析に携わっている方のご武運をお祈りいたします。

最後までお読みいただきありがとうございました。

参考