Decimal BASICで産業連関分析

十進BASIC/Decimal　BASICで産業連関分析をやってみる/やっている。

かれこれ20年になるがWindows時代になってF-BASICが有料になったのでこうなった。一言でいえばINV関数は便利、行列計算が1個ずつなのが難点、データの入出力もMAT文で読み込むときは行の最後にコンマをつけないといけない（らしい）。

私の勉強不足もあるかもしれないが、スカイライン図とか描画系となると、やはり強い。DBASIC（Decimal　BASICなのでこう略している）の強さは計算しながら図で確認できるところだ。
テキストの出力は遅く感じるけど、なにより言語の習得が圧倒的に楽。
講談社BlueBacksの木村良夫『パソコンを遊ぶ簡単プログラミング』１冊あれば習得できて数学も楽しめる。図が必須のシステム制御とかシステムダイナミックス的なこともできる。sampleコードの結構多いし、充実している。


下のコードは
内生344部門、粗付加価値10、需要50の395×395のデータを読み込んで逆行列を求め、横軸に算出額構成比、縦軸に産出乗数（レオンチェフ逆行列の列和）をとったグラフを描画するもの。

産業部門の規模と乗数を視覚化すれば、経済効果は高いが規模の小さい部門と、経済効果は低いが規模の大きな部門を視覚的に比較できるので、経済効果の誘発係数のみで議論する愚をさけることができる。

例えば、麦のように経済効果が高くても、経済規模が小さいような部門をどう評価するか。議論の出発点にはなるし、なぜ経済効果が高いのかという投入構造（費用構成）にまで議論を拡張できる。創始者のレオンチェフはなくなったが、産業連関分析の可能性はまだまだある。ちなみに三角化チャートの図も簡単にかけるし、ユニットストラクチャーも簡単に計算描画できる。取り出したい部門だけ描画するなど自由度も高いのがこの言語のいいところだ。



LET  n1=355
LET  n2=395
LET  n=344
LET  nn=50
LET  nnn=50
LET  mm=10
LET  mmm=10
dim datax(n1,n2)
DIM xi(n,n),v(mmm,n),f(n,nnn),m(n),x(n)
DIM a(n,n),mr(n),mc(n,n),mcv(n)
dim ii(n,n),xa(n,n),cb(n,n),sd(n),ob(n,n),ia(n,n)
rem **************** label ****************
DIM w$(n1),wv$(mmm),wfd$(n1)
OPEN #1:NAME "label_demand.txt"
FOR i=1 TO nnn
   INPUT #1:wfd$(i)
NEXT i
CLOSE #1
OPEN #1:NAME "label344.txt"
FOR i=1 TO n1
   INPUT #1:w$(i)
next i
CLOSE #1
FOR i=1 TO mmm
   LET wv$(i)=w$(n1-mmm+i)
NEXT i
REM **************** data（単位：百万円） ****************
LET f$="2005_okinawa_344.txt"
OPEN #1:NAME f$
MAT INPUT #1:datax
CLOSE #1
FOR i=1 TO n
   for j=1 to n
      LET  xi(i,j)=datax(i,j)
   next j
next i
REM MAT PRINT xi

for i=1 to n
   LET  x(i)=datax(n1,i)
NEXT i
REM MAT PRINT x

FOR i=1 TO n
   for j=1 to nnn
      LET  f(i,j)=datax(i,n+1+j)
   next j
next i
print "f"
REM MAT PRINT f

for i=1 to mmm
   for j=1 to n
      LET  v(i,j)=datax(n+1+i,j)
   next j
NEXT i

REM print "v"
rem FOR i=1 TO mmm
REM    PRINT wv$(i),
REM   FOR j=1 TO n
REM      PRINT v(i,j),
REM   NEXT j
REM    PRINT
rem NEXT i
REM MAT PRINT v
rem STOP

REM **************** 移輸入計383,384(344) ************
LET  importv=383
REM LET  transpv=96
FOR i=1 TO n
   REM PRINT (i),w$(i),datax(i,importv),datax(i,importv+1)
   LET  m(i)=datax(i,importv)+datax(i,importv+1)
next i
for k=1 to n
   LET  m(k)=-1*m(k)
next k
CLOSE #1
PRINT "m"
REM FOR i=1 TO n
REM   PRINT m(i)
REM next i
REM mat print m

for k=1 to n
   for kk=1 to n
      if x(kk)<>0 then LET  a(k,kk)=xi(k,kk)/x(kk) else LET  a(k,kk)=0
   next kk
next k

print "input coefficient a(i,j)"
REM mat print a
FOR i=1 TO n
   REM   PRINT w$(i),
   FOR j=1 TO n
      REM       PRINT a(i,j),
   NEXT j
   REM   PRINT
NEXT i

REM **************** 最終需要：地域内需要合計 ************
LET  fdemand=360
for k=1 to n
   LET  mr(k)=m(k)/(datax(k,fdemand))
   REM 移輸入額＞域内需要（定義的にあり得んと思うがとりあえず対策(^^;
   IF mr(k)>1 THEN LET mr(k)=1
   LET  mcv(k)=1-mr(k)
   REM PRINT k,mr(k),mcv(k),datax(k,fdemand),m(k),w$(k)
NEXT k
print "移輸入率（mr）"
REM mat print mr

FOR k=1 TO n
   for kk=1 to n
      if k=kk then LET  mc(k,kk)=1-mr(k)
      IF mc(k,kk)<0 k="" kk="" let="" mc="" p="" then="">   NEXT kk
NEXT k

print "自給率(I-M)"
REM print mc
PRINT "自給率（ベクトル）mcv"
rem MAT PRINT mcv

REM FOR i=1 TO n
REM    PRINT i,w$(i),mcv(i)
REM NEXT i

REM print "leontief inverse matrix"
mat ii=idn(n,n)
mat ia=ii-a
mat cb=inv(ia)
REM MAT PRINT cb

PRINT "出力の制御"

mat xa=mc*a
mat xa=ii-xa
REM mat print xa
mat ob=inv(xa)
REM print "inverse matrix  1/(I-(I-M)A)"
REM mat print ob

rem print "cb ob diagonal"
REM FOR i=1 TO n
REM    PRINT cb(i,i),ob(i,i),cb(i,i)-ob(i,i)
REM NEXT i

REM **********************************************************
REM **************  以上基本計算ルーチン  ********************
REM *******************    産出乗数   ************************
REM **********************************************************
REM 産出乗数、産出額構成比 x(i) -> cpx(i)
REM 横軸を需要や所得に変更して分析するのも可
REM
DIM multiple(n),cpx(n),multiplex(n)
FOR j=1 TO n
   LET summ=0
   LET summx=0
   FOR i=1 TO n
      LET  summ=summ+ob(i,j)
      LET  multiple(j)=summ
      LET  summx=summx+cb(i,j)
      LET  multiplex(j)=summx
   NEXT i
   PRINT j,multiple(j),multiplex(j),multiplex(j)-multiple(j),w$(j)
NEXT j
REM 構成比に用いる項目生産額:x(i)=datax(i,n1) ; 域内最終需要:datax(i,359),域内総需要(i,360),総生産(i,395)
REM 構成比は移輸出、移輸入、貿易の絶対額｛移輸出＋移輸入（絶対値）｝でも面白い

LET  s=0
FOR i=1 TO n
   REM  s=s+x(i)
   REM PRINT s
   LET  s=s+datax(i,395)
NEXT i
REM PRINT s
PRINT "部門名　輸入内生化型産出乗数　輸入外生化型算出乗数　構成比　面積"
FOR i=1 TO n
   REM  cpx(i)=100*x(i)/s
   LET  cpx(i)=100*datax(i,395)/s
   PRINT USING "#### " : i ;
   PRINT multiple(i),multiplex(i),cpx(i),100*multiple(i)*cpx(i),w$(i)
NEXT i

!'**********************************
dim sum(n)

LET  maxn=0
FOR i=1 TO n
   IF multiple(i)>=maxn THEN LET  maxn=multiple(i)
NEXT i

REM ******************************************
REM ************* 　部門数選択   *************
REM startsector、endsector：描画開始―終了部門
REM ******************************************
DIM dcpx(n),dmultiple(n),dmultiplex(n)
LET startsector=1
LET endsector=60
LET newn=endsector-startsector

FOR i=startsector TO endsector
   LET dcpx(i+1-startsector)=cpx(i)
   LET dmultiple(i+1-startsector)=multiple(i)
   LET dmultiplex(i+1-startsector)=multiplex(i)
NEXT i

LET n=newn
FOR i=1 TO n
   LET cpx(i)=dcpx(i)
   LET multiple(i)=dmultiple(i)
   LET multiplex(i)=dmultiplex(i)
NEXT i
LET s=0
FOR i=1 TO n
   LET s=s+dcpx(i)
NEXT i
FOR i=1 TO n
   LET  cpx(i)=100*dcpx(i)/s
NEXT i

PRINT "maxn : ",maxn
REM LET  maxn=4
REM ************** sect Window ***************
REM ******* SET WINDOW  left,right ,bottom,top
rem ************** 縦横軸補正 ****************
LET  botm=-1
LET  upper=6
REM *************** 右、左 *******************
LET  leftax=-2
LET  rightax=100
SET WINDOW leftax,rightax+2,botm-0.5,maxn+upper
PLOT LINES:leftax,0;rightax+2,0

PLOT TEXT ,AT 70,maxn+1:"Leontief Multiplier Chart : "
PLOT TEXT ,AT 70,maxn+2:f$

dim dx(n)
LET  s=0
for i=1 to n
   LET  s=s+cpx(i)
   LET  dx(i)=s
next i

dim pdot(n),edot(n)
LET  s=0
for i=1 to n
   LET  pdot(i)=dx(i)-(cpx(i)/1.6)
next i
dim numa(n),numb(n)
for i=1 to n step 2
   LET  numa(i)=i
   LET  numb(i)=i+1
   REM print numa(i),numb(i)
next i

LET  s=0
for i=1 to n
   LET  s=s+100/n
   LET  edot(i)=s-50/n
   IF i=numa(i) THEN PLOT TEXT ,AT edot(i)-0.5,0.8*botm:STR$(i+startsector) ELSE PLOT TEXT ,AT edot(i)-0.6,botm:STR$(i+startsector)
next i
SET LINE style 1
SET LINE WIDTH 1
SET LINE COLOR 1
for i=1 to n
   if i=numa(i) then plot lines:pdot(i),0;edot(i),0.5*botm else plot lines:pdot(i),0;edot(i),0.5*botm;edot(i),0.7*botm
next i

REM INPUT PROMPT "dum=":dummy
REM ************** revised ***************
REM prtxt Box内に表示する構成比の下限値、文字の表示位置 -> y:yst
SET LINE STYLE 1
SET LINE WIDTH 1
SET LINE COLOR 5
REM FOR i=1 TO n
REM PRINT i,cpx(i),dx(i),pdot(i),edot(i)
REM NEXT i
LET  prtxt=1.4
LET  yst=0.5
FOR i=1 TO n
   IF cpx(i)>=prtxt THEN  PLOT TEXT ,AT dx(i)-prtxt*cpx(i)/2, multiple(i)  :STR$(i+startsector)
NEXT i


rem *************** sector code ****************
REM 上の線：産出乗数の描画：輸入外生（標準乗数）
REM 必要が無ければこの線は描画しない

LET  sss=0
SET LINE STYLE 3
SET LINE WIDTH 1
SET LINE COLOR 4
plot 0,0;
FOR i=1 TO n
   PLOT LINES:sss,multiplex(i);
   LET  sss=sss+cpx(i)
   PLOT LINES:sss,multiplex(i);
   plot lines:sss,0;
   LET  sum(i)=sss
NEXT i

REM 上の線：産出乗数の描画：輸入内生（自給率調整乗数）
LET  sss=0
SET LINE STYLE 1
SET LINE WIDTH 1
SET LINE COLOR 1
plot 0,0;
for i=1 to n
   PLOT LINES:sss,multiple(i);
   LET  sss=sss+cpx(i)
   PLOT LINES:sss,multiple(i);
   plot lines:sss,0;
   LET  sum(i)=sss
NEXT i

REM ゼロの線：赤
PLOT LINES
SET LINE STYLE 1
SET LINE WIDTH 2
SET LINE COLOR 4
PLOT LINES:leftax,0;rightax+2,0

REM １の基線：青
PLOT LINES
SET LINE STYLE 3
SET LINE WIDTH 1
SET LINE COLOR 2
PLOT LINES:leftax,1;rightax+2,1



END