El Nino監視海域の海面温度データとSOIとの関係（クロス集計）
SOIと海面温度の関係
南方振動指数は大気海洋相互作用の「大気側」を見ている。
では海側はどうか？
何回も出てきた以下の図は、El Nino監視海域を示している。
（１）各監視海域の海面温度時系列とSOIとの相関関係
各海域の海面温度データ（DLすること）
NINO1,2,1+2,3,4,3.4の海面温度データでANOMは偏差を示す。
どの海域（どの列の海面温度データ）とSOIとの相関係数が大きくなるか調べよ。
（２）最小二乗法で求める回帰式
高い相関関係が得られた海域の海面温度を説明変数、SOIを説明変数として、単回帰式を求めよ。
Y（SOI）＝a × X(海面温度偏差） + b
の傾きaと切片bを求めるということ。
回帰式についてはここを参照
（３）決定係数
決定係数についてはここを参照
決定係数は回帰式がどれだけ当てはまっているかの指標。
以下はオリジナルのSOI時系列とSSTとの相関関係を調べるためのサンプルコード（ただし、いつも通り穴埋め形式）。
今回の課題は１２１か月移動平成分均を除去し、13か月移動平均をかけたSOI時系列とSSTを用いた相関係数、回帰式、決定係数の計算です。
c234567
character*120 header
real SOI(864),SOI2(502),date(864)
real YR(502),MO(502)
real nino(502,8),COR(8)
real eSOI(502),R2 !回帰式から求める水温、決定係数
!------- read SOI timeseries ---------!
inum=1
open(10,file="SOItimeseries.txt")
1000 read(10,*,end=555) date(inum),SOI(inum)
write(*,*) date(inum),SOI(inum)
inum=inum+1
goto 1000
555 inum=inum-1
close(10)
write(*,*) "number of data(SOI) = ",inum
!------- read SST timeseries ---------!
inum2=1
open(11,file="sstoi.indices.txt")
read(11,'(a120)') header
write(*,'(a120)') header
1001 read(11,*,end=666) YR(inum2),MO(inum2),(nino(inum2,k),k=1,8)
write(*,201) YR(inum2),MO(inum2),(nino(inum2,k),k=1,8)
inum2=inum2+1
goto 1001
666 inum2=inum2-1
close(11)
write(*,*) "number of data(SST) = ",inum2
201 format(f6.0,1x,f4.0,1x,8f8.2)
!------- find index to compare two datasets ---------!
! period: Jan. 1982 ~ Dec. 2022
!!２つのデータセットは期間が異なる。共通する期間の"座標"を見つける
!----------------------------------------------------!
istart=0 ; iend=0 ! SOI
istart2=0 ; iend2=0 ! SST
!!----SOI時系列で１９８２年１月と２０２２年１２月の場所を探す
icheck=1
do i=1,inum
if(icheck.eq.1.and.date(i).ge.1982)then
istart=i
icheck=icheck+1
endif
enddo
iend=inum
write(*,*) "SOI: istart,iend",istart,iend
write(*,*) "date start and end",date(istart),date(iend)
write(*,*) "icheck",icheck
!!----SST時系列で１９８２年１月と２０２２年１２月の場所を探す
icheck2=1
do j=1,inum2
if(変数YRが１９８２、かつ変数MOが１という条件)then
istart2=j
elseif(変数YRが２０２２、かつ変数MOが１２という条件)then
iend2=j
endif
enddo
write(*,*) "SST: istart2,iend2",istart2,iend2
write(*,*) "date start",YR(istart2),MO(istart2)
write(*,*) "date end",YR(iend2),MO(iend2)
!------------------------------------------------------
!------ Cal. correlation (SOI .vs. SST)
do j=1,8 ! choose region of SST monitoring area
! 1. get "MEAN" values !平均の計算
ave1=0 ! SOI
ave2=0 ! SST
idt=istart-istart2 !座標のギャップ
do i=istart,iend
ave1=ave1+SOI(i) ! sum of SOI
ave2=ave2+nino(i-idt,j)! sum of SST
enddo
dnum=iend-istart+1! データの数
ave1=ave1/dnum
ave2=ave2/dnum
!write(*,*) "idt,average",idt,ave1,ave2
! 2. get variance and std !分散・共分散と標準偏差
var1=0; var2=0;covar=0
do i=istart,iend ! SOIの座標を起点にしていることに注意
var1=var1+(SOI(i)-ave1)**2 ! sum of SOI anomaly
var2=var2+(nino(座標を計算,j)-ave2)**2 ! sum of SSTI anomaly
covar=covar+(SOI(i)-ave1)*(nino(座標を計算t,j)-ave2)
enddo
var1= !分散の計算
var2=!分散の計算
covar=!共分散の計算
std1=sqrt(var1)
std2=sqrt(var2)
! 3. cal correlation
!write(*,*) "var1,var2,covar",var1,var2,covar
!write(*,*) "std1,std2",std1,std2
COR(j)= !相関係数の計算
write(*,202) j,COR(j)
202 format("Area and Correlation-->",i4,f12.2)
enddo
!------ cal correlation finish! ----------------
!------------------------------------------------------
!======== LINEAR REGRESSION start! =============!
! X=SST(nino) & Y=SOI
!
do j=1,8 ! select sst data
B=0; C=0; D=0; E=0
do i=istart,iend
B=B+
C=C+
D=D+
E=E+
enddo
slope=
offset=
write(*,203) slope,offset
203 format("slope & offset -->",2f10.3)
C~~~~~~~~~~~~~~~~
C 決定係数の計算
C~~~~~~~~~~~~~~~~
!(1) SSTから予想されるSOIを回帰式から求める
do i=istart,iend
eSOI(i-idt)=
enddo
!(2) 残差平方和（SOIとeSOIの差を考える)
SSR=0
do i=istart,iend ! 373 ~ 864
SSR=SSR+
enddo
!(3) 全平方和TSSを求める
TSS=var1*dnum
R2=1-(SSR/TSS)
write(*,*) "SSR,TSS,R2",SSR,TSS,R2
enddo ! finish
stop
end