# バイオマス・エタノールのCO2削減効果分析
#
# バイマス・エタノールの投入係数，CO2原単位の計算は
# 岩井武士・熊谷真仁「バイオ起源のCO2排出と削減効果―ガソリン代替からダイエットまで―」
# 『三田商学研究学生論文集』2009年3月による
# 早見　均改定：104部門表への集計ほか
#
# 産業連関表は
# 総務省「平成12年産業連関表」104部門表（確報）2004年による
# http://www.stat.go.jp/data/io/index.htm
#
# CO2排出係数は
# 国立環境研究所地球環境研究センター
#「3EID産業連関表による環境負荷原単位データブック」2006年による。
# http://www-cger.nies.go.jp/publication/D031/jpn/page/data_file.htm
# 401部門表を104部門へ集計しCO2排出係数を使う
#
# データ              [単位]
# X,Fd,Ex,Im,VA,ct    [100万円]         104部門産業連関表
# bio$co2_g[1:104]    [t-CO2/100万円]   生産単位当たりCO2排出係数（ガソリン）
# bio$co2_ng[1:104]   [t-CO2/100万円]   生産単位当たりCO2排出係数（ガソリン以外）
# bio$co2_ng[105]     [t-CO2/t-Ethanol] バイオマス・エタノールの製造時のCO2排出係数
# bio$co2fd_g         [t-CO2/100万円]   ガソリン消費金額あたりガソリン燃焼時のCO2排出係数
# bio$co2fd_ng[1:104] [t-CO2/100万円]   ガソリン以外のエネルギー消費金額あたりのCO2排出係数
# bio$co2fd_ng[105]   [t-CO2/t-Ethanol] バイオマス・エタノール消費トンあたりCO2排出係数
# bio$Xgas_j          [100万円]         ガソリンの販路構成金額
# bio$Ai_ethanol      [t-CO2/t-Ethanol] バイオマス・エタノールの投入係数
# bio$fd_g            [100万円]         ガソリンの最終需要消費支出
# bio$fd_ng[1:104]    [100万円]         ガソリン以外の燃料の最終需要消費支出
# bio$fd_ng[105]      [t-Ethanol]       バイオマス・エタノールの消費重量トン
#

# データの読み込み
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)
bio<-read.table("BioEtOH.txt",header=TRUE,sep="\t")

# 計算用データの定義
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]
 }
#
# Simulation用データの定義
#
n_sim<-dim(bio)[1]
# 単位行列
Ident<-diag(1,nrow=n_sim,ncol=n_sim)
# 非消費支出最終需要 直接CO2が排出しない最終需要(政府消費には化石燃料支出はない)
f_nc<-vector(length=n_sim)
f_nc[1:n_sim]<-0
f_nc[1:n]<-gv+iv+ex-im
# ガソリン価格，対ガソリン・バイオマス燃費0.6，エタノール比重 0.77
gas_price<-87.892 # 2000年 Yen/Litre
i_gas<-29
i_etha<-105
gas_equiv<-1.0/0.6
etha_density<-0.77
# CO2排出係数 内生部門と最終需要 ガソリンとそれ以外
co2g<-as.vector(bio$co2_g)
co2n<-as.vector(bio$co2_ng)
co2fg<-as.vector(bio$co2fd_g)
co2fn<-as.vector(bio$co2fd_ng)
fg<-as.vector(bio$fd_g)
fn<-as.vector(bio$fd_ng)
#
# Simulation用関数の定義
#
get_A_ethanol<-function(A,X,ct,bio,gas_price,gas_equiv,etha_density,i_gas,i_etha,r_subst){
  n_bio<-dim(bio)[1]
  n_A<-dim(A)[1]
  A105 <- matrix(data=0,nrow=n_bio,ncol=n_bio)
  A105[1:n_A,1:n_A] <- A
  A105[1:n_A,i_etha] <- bio$Ai_ethanol[1:n_A]
  Xg<-as.vector(bio$Xgas_j)
  gas_x_litre<-Xg/gas_price # million litre
  etha_x_litre<-gas_x_litre*gas_equiv*r_subst
  for(j in 1:n_A){
    A105[i_gas,j]<-(X[i_gas,j]- Xg[j]*r_subst)/ct[j]
    A105[i_etha,j] <- etha_x_litre[j]*etha_density*1000.0/ct[j]
  }
  A105[i_etha,i_etha] <- 0
  return(A105)
}
get_fd_ethanol<-function(fg,fn,gas_price,gas_equiv,etha_density,i_gas,i_etha,r_subst){
  fd105<-fn
  fd105[i_gas]<-fd105[i_gas]+fg[i_gas]*(1-r_subst)
  gas_fd_litre<-(fg[i_gas]*r_subst)/gas_price # million litre
  etha_fd_litre<-gas_fd_litre*gas_equiv
  fd105[i_etha]<-etha_fd_litre*etha_density*1000.0 # tonne
  return(fd105)
}
#
# Simulation 00 : ethanol substitution (代替率) 0%
#
r_subst<-0.0
As00<-get_A_ethanol(A,X,x,bio,gas_price,gas_equiv,etha_density,i_gas,i_etha,r_subst)
fs00<-get_fd_ethanol(fg,fn,gas_price,gas_equiv,etha_density,i_gas,i_etha,r_subst)
Bs00<-solve(Ident - As00)
xs00<-Bs00%*%(f_nc+fs00)
CO2cn_00<-(co2n+co2g*(1-r_subst))*xs00
co2fd_s00<-co2fn*fn+co2fg*fg*(1-r_subst)
sum(CO2cn_00)+sum(co2fd_s00)
#
# Simulation 10 : ethanol substitution (代替率) 100%
#
r_subst<-1.00
As10<-get_A_ethanol(A,X,x,bio,gas_price,gas_equiv,etha_density,i_gas,i_etha,r_subst)
fs10<-get_fd_ethanol(fg,fn,gas_price,gas_equiv,etha_density,i_gas,i_etha,r_subst)
Bs10<-solve(Ident - As10)
xs10<-Bs00%*%(f_nc+fs10)
CO2cn_10<-(co2n+co2g*(1-r_subst))*xs10
co2fd_s10<-co2fg*fg*(1-r_subst)+co2fn*fn
sum(CO2cn_10)+sum(co2fd_s10)
#
# Simulation : ethanol substitution (代替率) 0-100%
#
j<-1
biz_cn<-numeric(101) # Carbon neutral
hhd_cn<-numeric(101) # Carbon neutral
biz_nn<-numeric(101) # adding CO2 from bioethanol
hhd_nn<-numeric(101) # adding CO2 from bioethanol
for( i in seq(1,101,by=5)){
  r_subst<-(i-1.0)/100.0
  As<-get_A_ethanol(A,X,x,bio,gas_price,gas_equiv,etha_density,i_gas,i_etha,r_subst)
  fs<-get_fd_ethanol(fg,fn,gas_price,gas_equiv,etha_density,i_gas,i_etha,r_subst)
  Bs<-solve(Ident - As)
  xs<-Bs%*%(f_nc+fs)
  biz_cn[j]<-sum((co2n+co2g*(1-r_subst))*xs)
  hhd_cn[j]<-sum(co2fn*fn+co2fg*fg*(1-r_subst))
  biz_nn[j]<-biz_cn[j]+co2fn[i_etha]*(xs[i_etha]-fs[i_etha])
  hhd_nn[j]<-hhd_cn[j]+sum(co2fn[i_etha]*fs[i_etha])
  j<-j+1
}
#
# グラフ表示用の操作
#
length(biz_cn)<-j-1
length(hhd_cn)<-j-1
length(biz_nn)<-j-1
length(hhd_nn)<-j-1

dev.off()
xlbl<-seq(0,101,by=5)
y_max<-max(biz_cn+hhd_cn)
plot(biz_cn,xlab="代替率 (%)",ylab="総CO2排出量 (百万t-CO2)",axes=0,type="l",ylim=c(1100*10^6,y_max))
lines(biz_cn+hhd_cn)
axis(1, 1:21,labels = formatC(xlbl, format="fg"))
(my.at<-seq(10^9*0.95,10^9*1.36,,by=10^7*5))
axis(2, my.at,labels = formatC(my.at/10^6, format="fg"))
text(15,1120*10^6,"(産業部門)")
text(5,1200*10^6,"(家計部門)")
text(15,1300*10^6,"国内排出量")
cbind(biz_cn,hhd_cn,biz_cn+hhd_cn)

dev.new() # 別のWindowsに表示する場合
y_max<-max(biz_nn+hhd_nn)
plot(biz_nn,xlab="代替率 (%)",ylab="総CO2排出量 (百万t-CO2)",axes=0,type="l",ylim=c(1100*10^6,y_max),main="バイオマス・エタノールからのCO2を加えた場合")
lines(biz_nn+hhd_nn)
lines(biz_cn,col="Red")
lines(biz_cn+hhd_cn,col="Red")
axis(1, 1:21,labels = formatC(xlbl, format="fg"))
(my.at<-seq(10^9*0.95,10^9*1.46,,by=10^7*5))
axis(2, my.at,labels = formatC(my.at/10^6, format="fg"))
text(15,1120*10^6,"(産業部門)")
text(5,1250*10^6,"(家計部門)")
text(15,1325*10^6,"国内排出量")
