# 産業連関分析・競争輸入型
# 2000年: 総務省「平成12年産業連関表」104部門表による
# 104部門へ集計されたCO2排出係数を使う
#
# データ          [単位]
# co2_coef        [CO2-ton/100万円]
# co2_104         [CO2-ton]
# ct_104          [100万円]
# co2_nh_direct <-1395770   # 家計外消費の燃料消費からのCO2
# co2_hh_heating<-92663224  # 家計の暖房厨房(都市ガス・石炭・灯油・LPG)からのCO2
# co2_hh_trans  <-82900553  # 家計の自動車(ガソリン・軽油)からのCO2
#

# データの読み込み
load("co2.Rdata")
X<-read.table("X104.txt",header=TRUE,sep="\t")
Fd<-read.table("fd104.txt",header=TRUE,sep="\t")
Ex<-read.table("ex104.txt",header=TRUE,sep="\t")
Im<-read.table("im104.txt",header=TRUE,sep="\t")
VA<-read.table("v104.txt",header=TRUE,sep="\t")
ct<-read.table("ct104.txt",header=TRUE)

# 計算用データの定義
n<-dim(X)[1]
I<-diag(1,nrow=n,ncol=n)
A<-array(data=0,dim=c(n,n))
cs<-array(data=0,dim=c(n))
gv<-array(data=0,dim=c(n))
iv<-array(data=0,dim=c(n))
ex<-array(data=0,dim=c(n))
x<-array(data=0,dim=c(n))
va<-array(data=0,dim=c(n))
v<-array(data=0,dim=c(n))

# 家計外消費支出+民間消費支出
cs<-Fd[,1]+Fd[,2]
# 一般政府最終消費支出+固定資本形成(公的)
gv<-Fd[,3]+Fd[,4]+Fd[,5]
# 固定資本形成(民間)+在庫純増
iv<-Fd[,6]+Fd[,7]
# 輸出+調整項
ex<-Ex[,1]+Ex[,2]
# 輸入+関税+輸入商品税: 輸入データはもともとマイナスで入力されている．
im<--(Im[,1]+Im[,2]+Im[,3])
# 粗付加価値計
va<-VA[1,]+VA[2,]+VA[3,]+VA[4,]+VA[5,]+VA[6,]+VA[7,]
# 国内生産額
x<-ct[,1]

# 投入係数行列の計算
for( j in c(1:n)) {
  A[,j]<-X[,j]/x[j]
 }

# 輸入係数行列の計算
imm<-im/(x+im-ex)
M<-diag(imm)

# レオンチェフ逆行列の計算
BM<-solve(I-(I-M)%*%A)

F<-as.matrix(cbind(Fd[,1:5],iv))
colnames(F)<-c(colnames(Fd[,1:5]),"国内総固定資本形成.民間")

# 生産誘発額 : Xfd
Xfd<-BM%*%as.matrix(cbind((I-M)%*%F,ex))
colnames(Xfd)<-c(colnames(F),"輸出")

# 最終需要で誘発されるCO2排出量
CO2fd<-co2_coef%*%Xfd
CO2fd_vol<-sprintf("%.2f億t",CO2fd/10^8)
lbl<-sprintf("%s%s",colnames(CO2fd),CO2fd_vol)

# CO2 induced by final demand
pie(CO2fd,labels=lbl,main="最終需要で誘発されるCO2排出量")

dev.new()
CO2fd_d<-cbind(CO2fd,co2_nh_direct,co2_hh_heating,co2_hh_trans)
colnames(CO2fd_d)<-c(colnames(CO2fd),"家計外消費(直)","家計暖房・厨房(直)","家計自動車(直)")
CO2fd_v_d<-sprintf("%.2f億t",CO2fd_d/10^8)
lbl<-sprintf("%s%s",colnames(CO2fd_d),CO2fd_v_d)
pie(CO2fd_d,labels=lbl,main="最終需要の直接排出と誘発されたCO2排出量")

# 炭素税の転嫁
ell<-c(1:n)
ell[1:n]<-1
tax_c<-t(BM)%*%diag(co2_coef)%*%ell
j_sec<-match(sort(tax_c,decreasing=TRUE),tax_c)
cbind(j_sec[1:20],tax_c[j_sec[1:20]])
