github.com/Anderson-Lu/gobox@v0.0.0-20191127065433-3e6c4c2da420/math/adf/adf.go (about)

     1  package adf
     2  
     3  import (
     4  	"github.com/berkmancenter/ridge"
     5  	"github.com/gonum/matrix/mat64"
     6  	"github.com/gonum/stat"
     7  
     8  	"math" 
     9  )
    10  
    11  const (
    12  	LPenalty      = 0.0001 // L penalty to pass to ridge regression
    13  	DefaultPValue = -3.45  // Test p-value threshold
    14  )
    15  
    16  // An instance of an ADF test
    17  type ADF struct {
    18  	Series          []float64 // The time series to test
    19  	PValueThreshold float64   // The p-value threshold for the test
    20  	Statistic       float64   // The test statistic
    21  	Lag             int       // The lag to use when running the test
    22  }
    23  
    24  // New creates and returns a new ADF test.
    25  func New(series []float64, pvalue float64, lag int) *ADF {
    26  	if pvalue == 0 {
    27  		pvalue = DefaultPValue
    28  	}
    29  	if lag < 0 {
    30  		lag = int(math.Floor(math.Cbrt(float64(len(series)))))
    31  	}
    32  	newSeries := make([]float64, len(series))
    33  	copy(newSeries, series)
    34  	return &ADF{Series: newSeries, PValueThreshold: pvalue, Lag: lag}
    35  }
    36  
    37  // Run runs the Augmented Dickey-Fuller test.
    38  func (adf *ADF) Run() {
    39  	series := adf.Series
    40  	mean := stat.Mean(series, nil)
    41  	if mean != 0.0 {
    42  		for i, v := range series {
    43  			series[i] = v - mean
    44  		}
    45  	}
    46  	n := len(series) - 1
    47  	y := diff(series)
    48  	lag := adf.Lag
    49  	k := lag + 1
    50  	z := laggedMatrix(y, k)
    51  	zcol1 := mat64.Col(nil, 0, z)
    52  	xt1 := series[k-1 : n]
    53  	r, c := z.Dims()
    54  	var design *mat64.Dense
    55  
    56  	if k > 1 {
    57  		yt1 := z.View(0, 1, r, c-1)
    58  		design = mat64.NewDense(n-k+1, k, nil)
    59  		design.SetCol(0, xt1)
    60  
    61  		_, c = yt1.Dims()
    62  		for i := 0; i < c; i++ {
    63  			design.SetCol(1+i, mat64.Col(nil, i, yt1))
    64  		}
    65  	} else {
    66  		design = mat64.NewDense(n-k+1, 1, nil)
    67  		design.SetCol(0, xt1)
    68  	}
    69  	regressY := mat64.NewVector(len(zcol1), zcol1)
    70  
    71  	rr := ridge.New(design, regressY, LPenalty)
    72  	rr.Regress()
    73  	beta := rr.Coefficients.RawVector().Data
    74  	sd := rr.StdErrs
    75  
    76  	adf.Statistic = beta[0] / sd[0]
    77  }
    78  
    79  // IsStationary returns true if the tested time series is stationary.
    80  func (adf ADF) IsStationary() bool {
    81  	return adf.Statistic < adf.PValueThreshold
    82  }
    83  
    84  func diff(x []float64) []float64 {
    85  	y := make([]float64, len(x)-1)
    86  	for i := 0; i < len(x)-1; i++ {
    87  		y[i] = x[i+1] - x[i]
    88  	}
    89  	return y
    90  }
    91  
    92  func laggedMatrix(series []float64, lag int) *mat64.Dense {
    93  	r, c := len(series)-lag+1, lag
    94  	m := mat64.NewDense(r, c, nil)
    95  	for j := 0; j < c; j++ {
    96  		for i := 0; i < r; i++ {
    97  			m.Set(i, j, series[lag-j-1+i])
    98  		}
    99  	}
   100  	return m
   101  }