scalaで単位根検定 unit root test (Phillips-Perron Tests in scala)

Scalaで単位根検定　unit root test (Phillips-Perron Tests in scala)
RのPP.testをscalaに翻訳

import scala.math._
import org.apache.commons.math3.stat.correlation.SpearmansCorrelation;
import org.apache.commons.math3.stat.regression.OLSMultipleLinearRegression;
import org.apache.commons.math3.linear
import org.apache.commons.math3.analysis.interpolation.SplineInterpolator
import org.apache.commons.math3.analysis.interpolation.LinearInterpolator
object PP_test{


        def PP_test_func(x:Array[Double]):Double={
                var y=1.0;
                if( x.length <= 1){ println("x is not a vector");return -1.0;}

                var yt=x.slice(1,x.length)

                var yt1=x.slice(0,(x.length-1))
                var n=yt.length
                var tt=(1 to n).map{a=> a-n/2.0}.toArray
                var reg=new OLSMultipleLinearRegression();

                var y1=(yt)

                val x1=Array.ofDim[Double](n,2)
                (0 to (n-1)).foreach{a=>
                        x1(a)(0)=tt(a);
                        x1(a)(1)=yt1(a);


                }


                reg.newSampleData(y1,x1);
                var beta=(reg.estimateRegressionParameters());
                var beta_v=reg.estimateRegressionParametersVariance();
                var epsilon=reg.estimateRegressionStandardError();

                var residuals=reg.estimateResiduals();
                println("test"+beta(0)+" "+beta(1)+" "+beta(2)+" "+sqrt(epsilon*epsilon*beta_v(2)(2)));             println("var");


                println("epsilon");

                residuals.foreach{println}

                //Need Try and catch

                var tmp_int=0;
                var tstat=(beta(2)-1)/sqrt(epsilon*epsilon*beta_v(2)(2));          var  ssqru=residuals.map{a=> a*a}.sum/n;

                residuals.foreach{println}
                var l=(4.0*sqrt(sqrt(n/100.0))).toInt
                tmp_int=0;
                var gamma2=(1 to (l)).map{j=> var k=0.0; ((j+1) to (n)).foreach{t=> k=k+residuals(t-1)*residuals(t-j-1);};k/n; } ;
                println("gamma2 ", gamma2(0));
                tmp_int=0;
                var ssqrtl=ssqru+2.0*gamma2.map{a=> tmp_int=tmp_int+1;(1-(tmp_int)/(1.0+l))*a}.sum;
                var n2=n*n;

                var trm1=n2*(n2-1.0)*(yt1.map{a=> a*a}.sum)/12.0;
                tmp_int=0;
                var trm2=n*pow(yt1.map{a=> tmp_int=tmp_int+1;a*tmp_int}.sum,2);
                tmp_int=0;
                var trm3=n*(n+1.0)*yt1.map{a=> tmp_int=tmp_int+1;a*tmp_int}.sum*yt1.sum;

                var trm4=(n*(n+1.0)*(2.0*n+1)*pow(yt1.sum,2))/6.0
                var Dx=trm1-trm2+trm3-trm4;
                var STAT=sqrt(ssqru)/sqrt(ssqrtl)*tstat-pow(n,3)/(4.0*sqrt(3)*sqrt(Dx)*sqrt(ssqrtl))*(ssqrtl-ssqru);


                println(trm1," ",trm2," ",trm3," ",trm4," ",Dx," ",STAT," ",ssqrtl);

                var table=Array(Array(4.38, 4.15, 4.04, 3.99, 3.98, 3.96), Array(3.95, 3.8, 3.73, 3.69, 3.68, 3.66), Array(3.6, 3.5, 3.45, 3.43, 3.42, 3.41), Array(3.24, 3.18, 3.15, 3.13, 3.13, 3.12), Array(1.14, 1.19, 1.22, 1.23, 1.24, 1.25), Array(0.8, 0.87, 0.9, 0.92, 0.93, 0.94), Array(0.5, 0.58, 0.62, 0.64, 0.65, 0.66), Array(0.15, 0.24, 0.28, 0.31, 0.32, 0.33))

                table=table.map{a=> a.map{-_}}
                //table.foreach{a=> a.foreach{println}}
                var tablen=table.length;
                var tableT=Array(25.0,50.0,100.0,250,500,100000)
                var tablep=Array(0.01,0.025,0.05,0.1,0.9,0.95,0.975,0.99)
                //var tableipl: Array[Double]=new Array[Double](n);
                var tableipl=Array.fill(n)(0.0);

                //tableipl.foreach{println}


                //var interpolator=new SplineInterpolator();

                var interpolator=new LinearInterpolator();
                var m=0;

                tableipl=table.map{a:Array[Double] =>

                        var func=interpolator.interpolate(tableT,a);
                        if(n<25) m=25 else m=n;
                        if(m>100000) m=100000 else m=m;
                        func.value(m+0.0);

                }
                tableipl.foreach{println}

                var STAT2=0.0
                if(STAT< tableipl.min) STAT2=tableipl.min else STAT2=STAT;
                if(STAT> tableipl.max) STAT2=tableipl.max else STAT2=STAT2;

                var PVAL=interpolator.interpolate(tableipl,tablep).value(STAT2);

                println("PVAL:",PVAL);


                PVAL;
        }



        def main(args: Array[String]){
                val data1=Array(20,17,30,42.3,17,50);
                

                println(PP_test_func(data1));

        }
}

ファイルからデータを読み込む

import scala.math._
import org.apache.commons.math3.stat.correlation.SpearmansCorrelation;
import org.apache.commons.math3.stat.regression.OLSMultipleLinearRegression;
import org.apache.commons.math3.linear
import org.apache.commons.math3.analysis.interpolation.SplineInterpolator
import org.apache.commons.math3.analysis.interpolation.LinearInterpolator
import scala.io.Source
object PP_test2{



        def PP_test_func(x:Array[Double]):Double={
                var y=1.0;
                if( x.length <= 1){ println("x is not a vector");return -1.0;}

                var yt=x.slice(1,x.length)

                var yt1=x.slice(0,(x.length-1))
                var n=yt.length
                var tt=(1 to n).map{a=> a-n/2.0}.toArray
                var reg=new OLSMultipleLinearRegression();

                var y1=(yt)

                val x1=Array.ofDim[Double](n,2)
                (0 to (n-1)).foreach{a=>
                        x1(a)(0)=tt(a);
                        x1(a)(1)=yt1(a);


                }


                reg.newSampleData(y1,x1);
                var beta=(reg.estimateRegressionParameters());
                var beta_v=reg.estimateRegressionParametersVariance();
                var epsilon=reg.estimateRegressionStandardError();

                var residuals=reg.estimateResiduals();
                println("test"+beta(0)+" "+beta(1)+" "+beta(2)+" "+sqrt(epsilon*epsilon*beta_v(2)(2)));             println("var");


                println("epsilon");

                residuals.foreach{println}

                //Need Try and catch

                var tmp_int=0;
                var tstat=(beta(2)-1)/sqrt(epsilon*epsilon*beta_v(2)(2));          var  ssqru=residuals.map{a=> a*a}.sum/n;

                residuals.foreach{println}
                var l=(4.0*sqrt(sqrt(n/100.0))).toInt
                tmp_int=0;
                var gamma2=(1 to (l)).map{j=> var k=0.0; ((j+1) to (n)).foreach{t=> k=k+residuals(t-1)*residuals(t-j-1);};k/n; } ;
                println("gamma2 ", gamma2(0));
                tmp_int=0;
                var ssqrtl=ssqru+2.0*gamma2.map{a=> tmp_int=tmp_int+1;(1-(tmp_int)/(1.0+l))*a}.sum;
                var n2=n*n;

                var trm1=n2*(n2-1.0)*(yt1.map{a=> a*a}.sum)/12.0;
                tmp_int=0;
                var trm2=n*pow(yt1.map{a=> tmp_int=tmp_int+1;a*tmp_int}.sum,2);
                tmp_int=0;
                var trm3=n*(n+1.0)*yt1.map{a=> tmp_int=tmp_int+1;a*tmp_int}.sum*yt1.sum;

                var trm4=(n*(n+1.0)*(2.0*n+1)*pow(yt1.sum,2))/6.0
                var Dx=trm1-trm2+trm3-trm4;
                var STAT=sqrt(ssqru)/sqrt(ssqrtl)*tstat-pow(n,3)/(4.0*sqrt(3)*sqrt(Dx)*sqrt(ssqrtl))*(ssqrtl-ssqru);


                println(trm1," ",trm2," ",trm3," ",trm4," ",Dx," ",STAT," ",ssqrtl);

                var table=Array(Array(4.38, 4.15, 4.04, 3.99, 3.98, 3.96), Array(3.95, 3.8, 3.73, 3.69, 3.68, 3.66), Array(3.6, 3.5, 3.45, 3.43, 3.42, 3.41), Array(3.24, 3.18, 3.15, 3.13, 3.13, 3.12), Array(1.14, 1.19, 1.22, 1.23, 1.24, 1.25), Array(0.8, 0.87, 0.9, 0.92, 0.93, 0.94), Array(0.5, 0.58, 0.62, 0.64, 0.65, 0.66), Array(0.15, 0.24, 0.28, 0.31, 0.32, 0.33))

                table=table.map{a=> a.map{-_}}
                //table.foreach{a=> a.foreach{println}}
                var tablen=table.length;
                var tableT=Array(25.0,50.0,100.0,250,500,100000)
                var tablep=Array(0.01,0.025,0.05,0.1,0.9,0.95,0.975,0.99)
                //var tableipl: Array[Double]=new Array[Double](n);
                var tableipl=Array.fill(n)(0.0);

                //tableipl.foreach{println}


                //var interpolator=new SplineInterpolator();

                var interpolator=new LinearInterpolator();
                var m=0;

                tableipl=table.map{a:Array[Double] =>

                        var func=interpolator.interpolate(tableT,a);
                        if(n<25) m=25 else m=n;
                        if(m>100000) m=100000 else m=m;
                        func.value(m+0.0);

                }
                tableipl.foreach{println}

                var STAT2=0.0
                if(STAT< tableipl.min) STAT2=tableipl.min else STAT2=STAT;
                if(STAT> tableipl.max) STAT2=tableipl.max else STAT2=STAT2;

                var PVAL=interpolator.interpolate(tableipl,tablep).value(STAT2);

                println("PVAL:",PVAL);


                PVAL;
        }



        def main(args: Array[String]){
              
                var source1 =Source.fromFile("test.csv","utf-8")


                val lines1=source1.getLines
                var list1=List[Double]()
                lines1.foreach{a:String => list1=(a.split('\t')(4).toDouble)::list1;}
                var rev_list1=list1.reverse;
               
                println(PP_test_func(rev_list1.toArray));

        }
}

そのほか 移動法など

import scala.io.Source
import org.apache.commons.math3.stat.correlation.SpearmansCorrelation

import scala.math._
import org.apache.commons.math3.stat.correlation.SpearmansCorrelation;
import org.apache.commons.math3.stat.regression.OLSMultipleLinearRegression;
import org.apache.commons.math3.linear
import org.apache.commons.math3.analysis.interpolation.SplineInterpolator
import org.apache.commons.math3.analysis.interpolation.LinearInterpolator
object cor_cehek2_scala{



        def PP_test(x:Array[Double]):Double={
                var y=1.0;
                if( x.length <= 1){ println("x is not a vector");return -1.0;}

                var yt=x.slice(1,x.length)

                var yt1=x.slice(0,(x.length-1))
                var n=yt.length
                var tt=(1 to n).map{a=> a-n/2.0}.toArray
                var reg=new OLSMultipleLinearRegression();

                var y1=(yt)

                val x1=Array.ofDim[Double](n,2)
                (0 to (n-1)).foreach{a=>
                        x1(a)(0)=tt(a);
                        x1(a)(1)=yt1(a);


                }


                reg.newSampleData(y1,x1);
                var beta=(reg.estimateRegressionParameters());
                var beta_v=reg.estimateRegressionParametersVariance();
                var epsilon=reg.estimateRegressionStandardError();

                var residuals=reg.estimateResiduals();
                //println("test"+beta(0)+" "+beta(1)+" "+beta(2)+" "+sqrt(epsilon*epsilon*beta_v(2)(2)));           println("var");


                //println("epsilon");

                //residuals.foreach{println}

                //Need Try and catch

                var tmp_int=0;
                var tstat=(beta(2)-1)/sqrt(epsilon*epsilon*beta_v(2)(2));                   var ssqru=residuals.map{a=> a*a}.sum/n;

                //residuals.foreach{println}
                var l=(4.0*sqrt(sqrt(n/100.0))).toInt
                tmp_int=0;
                var gamma2=(1 to (l)).map{j=> var k=0.0; ((j+1) to (n)).foreach{t=> k=k+residuals(t-1)*residuals(t-j-1);};k/n; } ;
                //println("gamma2 ", gamma2(0));
                tmp_int=0;
                var ssqrtl=ssqru+2.0*gamma2.map{a=> tmp_int=tmp_int+1;(1-(tmp_int)/(1.0+l))*a}.sum;
                var n2=n*n;

                var trm1=n2*(n2-1.0)*(yt1.map{a=> a*a}.sum)/12.0;
                tmp_int=0;
                var trm2=n*pow(yt1.map{a=> tmp_int=tmp_int+1;a*tmp_int}.sum,2);
                tmp_int=0;
                var trm3=n*(n+1.0)*yt1.map{a=> tmp_int=tmp_int+1;a*tmp_int}.sum*yt1.sum;

                var trm4=(n*(n+1.0)*(2.0*n+1)*pow(yt1.sum,2))/6.0
                var Dx=trm1-trm2+trm3-trm4;
                var STAT=sqrt(ssqru)/sqrt(ssqrtl)*tstat-pow(n,3)/(4.0*sqrt(3)*sqrt(Dx)*sqrt(ssqrtl))*(ssqrtl-ssqru);


                //println(trm1," ",trm2," ",trm3," ",trm4," ",Dx," ",STAT," ",ssqrtl);

                var table=Array(Array(4.38, 4.15, 4.04, 3.99, 3.98, 3.96), Array(3.95, 3.8, 3.73, 3.69, 3.68, 3.66), Array(3.6, 3.5, 3.45, 3.43, 3.42, 3.41), Array(3.24, 3.18, 3.15, 3.13, 3.13, 3.12), Array(1.14, 1.19, 1.22, 1.23, 1.24, 1.25), Array(0.8, 0.87, 0.9, 0.92, 0.93, 0.94), Array(0.5, 0.58, 0.62, 0.64, 0.65, 0.66), Array(0.15, 0.24, 0.28, 0.31, 0.32, 0.33))

                table=table.map{a=> a.map{-_}}
                //table.foreach{a=> a.foreach{println}}
                var tablen=table.length;
                var tableT=Array(25.0,50.0,100.0,250,500,100000)
                var tablep=Array(0.01,0.025,0.05,0.1,0.9,0.95,0.975,0.99)
                //var tableipl: Array[Double]=new Array[Double](n);
                var tableipl=Array.fill(n)(0.0);

                //tableipl.foreach{println}


                //var interpolator=new SplineInterpolator();

                var interpolator=new LinearInterpolator();
                var m=0;

                tableipl=table.map{a:Array[Double] =>

                        var func=interpolator.interpolate(tableT,a);
                        if(n<25) m=25 else m=n;
                        if(m>100000) m=100000 else m=m;
                        func.value(m+0.0);

                }
                //tableipl.foreach{println}

                var STAT2=0.0
                if(STAT< tableipl.min) STAT2=tableipl.min else STAT2=STAT;
                if(STAT> tableipl.max) STAT2=tableipl.max else STAT2=STAT2;

                var PVAL=interpolator.interpolate(tableipl,tablep).value(STAT2);

                //println("PVAL:",PVAL);


                PVAL;
        }





         def box_mean(x:List[Double],n:Int=2):List[Double]={

                var y=List[Double]();
                var k:Int=0;
                var nn=n-1;
                var sum=0.0;

                x.foreach{ a =>
                        k match{
                                case k if k%n==n-1 => sum=sum+a/n;y=sum::y;sum=0;
                                case _ => sum=sum+a/n;
                        }

                        k=k+1;
                }
               y.reverse
        }



         def diff_boxmean(x:List[Double],n:Int=2):List[Double]={
                var tmp=0.0;
                var y=List[Double]();
                box_mean(x,n).foreach{a =>

                        y=a-tmp::y;
                        tmp=a;

                }
                y.reverse.tail
        }

        def diff_boxmean_cor(x1:List[Double],x2:List[Double],n:Int=2):Double={
                var xx1=diff_boxmean(x1,n)
                var xx2=diff_boxmean(x2,n)
                val cor = new SpearmansCorrelation();
                cor.correlation(xx1.toArray,xx2.toArray);
        }



        def diff_boxmean_cor_PPtest(x1:List[Double],x2:List[Double],n:Int=2):Array[Double]={
                var xx1=diff_boxmean(x1,n)
                var xx2=diff_boxmean(x2,n)
                val cor = new SpearmansCorrelation();
                var cor2=cor.correlation(xx1.toArray,xx2.toArray);
                var pp1=PP_test(xx1.toArray)

                var pp2=PP_test(xx2.toArray)

                Array(cor2,pp1,pp2)
        }




        def main(args: Array[String]){




                var source1 =Source.fromFile("data1.txt","utf-8")
                var source2 =Source.fromFile("data2.txt","utf-8")
                var source3 =Source.fromFile("normalize2.txt","utf-8")
                var source4 =Source.fromFile("normalize1.txt","utf-8")
     
                val lines1=source1.getLines
                val lines2=source2.getLines
                val lines3=source3.getLines
                val lines4=source4.getLines


                var list1=List[Double]();
                var list2=List[Double]();
                var list3=List[Double]();
                var list4=List[Double]();
              
                lines1.foreach{a:String => list1=(a.split('\t')(4).toDouble)::list1;}
                lines2.foreach{a:String => list2=(a.split(' ')(1).toDouble)::list2;}

                lines3.foreach{a:String => list3=(a.toDouble)::list3;}

                lines4.foreach{a:String => list4=(a.toDouble)::list4;}
                var rev_list1=list1.reverse;
                list2=list2.tail.tail.tail;
                var rev_list2=list2.reverse;
                var rev_list3=list3.reverse;
                var rev_list4=list4.reverse;

            
                rev_list2=rev_list2.map{a =>
                        var c= a/rev_list3.head
                        rev_list3=rev_list3.tail;
                        c
                }

                rev_list1=rev_list1.map{a =>
                        var c= a/rev_list4.head
                        rev_list4=rev_list4.tail;
                        c
                }

  
                source1.close();
                source2.close();



                //exit();
                var cor1=new Array[Double](150);
                var p1=new Array[Double](150);
                var p2=new Array[Double](150);
                var pp=new Array[Double](3);
                for(i <-1 to 150){
                        println(i+" "+diff_boxmean_cor_PPtest(rev_list1,rev_list2,i).mkString(" "));

                        pp=diff_boxmean_cor_PPtest(rev_list1,rev_list2,i)
                        cor1(i-1)=pp(0)
                        p1(i-1)=pp(1)
                        p2(i-1)=pp(2)

                }


                var minmin= Array(p1.indexWhere(a=> a>=0.05) match {case -1 =>1000; case other =>other} ,p2.indexWhere(a=> a>=0.05)).min
                println(minmin);
       
        }
}