Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm trying to create a function that calculates and plots the sample variogram for a time series. It should take two arguments: data, a time series, essentially a list of reals; k, an integer greater than 0. It should produce a ListLinePlot with $k$ on the horizontal axis and $\hat{G_k}$ on the vertical axis.

The plotting and such isn't really the problem. The problem is getting the mathematical formulas, which are below, into Mathematica formulas.

$\hat{G_k}$ is essentially a ratio of variance at lag $k$ over the variance at lag $1$:

$\hat{G_k}=\frac{s^2_k}{s^2_1},k=1,2,...$

$s^2_k=\frac{\sum_{t=1}^{n-k}(d^k_t-\bar{d}^k)^2}{n-k-1}$

$\bar{d}^k=(n-k)^{-1}\sum{}d^k_t$

$d^k_t=z_{t+k}-z_t$

$z_t$ is the time series of the length $n$.

I've been beating my head against this for a while now. Any help would be appreciated.

Edit: A variogram is statistical tool to assess whether a time series is stationary or not. For a stationary it starts out at $1$ for $k=1$ and then increases until it starts to hover around a certain value. The theory and formulas I have here is from Time Series Analysis and Forecasting by Example. You might be able to use the Search Inside at Amazon to have a look at it.

Edit2: Here is a demoplot based on dummydata as to how the end result should look for a stationary time series. Mathematica graphics

share|improve this question
1  
What's a variogram? Can you provide a link? –  rcollyer Feb 28 '12 at 19:02
    
rcollyer: It's a version of a correlogram and can be directly computed with ListConvolve, ListCorrelate, et al. –  whuber Feb 28 '12 at 19:39
    
As ListPlot can do the same thing as ListLinePlot, I removed the version-8 tag. –  rcollyer Feb 28 '12 at 19:42
    
@rcollyer ListLinePlot was "New in 6" -- not a problem even for me ;-) –  Mr.Wizard Feb 28 '12 at 20:38
    
@Mr.Wizard I swore it was new in 8. –  rcollyer Feb 28 '12 at 22:33

4 Answers 4

up vote 2 down vote accepted

Lag differences are not the same as second-differencing etc, so Differences is not the right approach.

test2 = Array[f, 10]

In[23]:= Differences[test2, 3]

Out[23]= {-f[1] + 3 f[2] - 3 f[3] + f[4], -f[2] + 3 f[3] - 3 f[4] + f[5], -f[3] + 3 f[4] - 3 f[5] + f[6], -f[4] + 3 f[5] - 3 f[6] + f[7], -f[5] + 3 f[6] - 3 f[7] + f[8], -f[6] + 3 f[7] - 3 f[8] + f[9], -f[7] + 3 f[8] - 3 f[9] + f[10]}

Building on kguler's answer, here is an alternative way of getting lag-differences. I have written it up as a separate function because it has more general utility than this specific application.

lagdif[l_List, k_Integer?Positive] /; k < Length[l] := 
  Drop[l, k] - Drop[l, -k]

We can then use a similar approach:

dif[list_] := Table[lagdif[list, k], {k, Length[list] - 1}];
dbar[list_] := Mean /@ dif[list];
var[list_] := Variance /@ (Most@(dif[list]));
variogram[list_] := ((Rest@#)/First@#) &@var[list];

Testing with I(1) data (i.e. non-stationary with a single unit root):

$y_t = y_{t-1} + \epsilon_t$

data = Accumulate[RandomReal[{-1, 1}, 20]];

GraphicsRow[{ListLinePlot[data], ListLinePlot[variogram[data]]}]

enter image description here

This makes me wonder how powerful the variogram is for data with shocks that have bounded support.

share|improve this answer

I don't know if it is OK to answer your own question, but I thought I would give an answer as to what I ultimately ended up with. This is essentially based on Verbeia's answer.

I've added a option for how many lags to calculate since calculating the variance for lags up to $n-1$ doesn't make much sense, since the number of data points you have decrease with the size of the lag. This leads to the estimate of the variance being uncertain.

lagdif[l_List, k_Integer?Positive] /; k < Length[l] := 
 Drop[l, k] - Drop[l, -k]
variogram[list_, 
   k_] := (#/First@#) &@(Variance /@ Table[lagdif[list, i], {i, k}]);
plotvariogram[list_, k_] := 
 GraphicsRow[{ListPlot[list, Joined -> True, 
    AxesLabel -> {"t", "\!\(\*SubscriptBox[\(X\), \(t\)]\)"}], 
   ListPlot[variogram[list, k], AxesOrigin -> {0, 0}, 
    AxesLabel -> {"k", 
      "\!\(\*SubscriptBox[OverscriptBox[\(G\), \(^\)], \(k\)]\)"}, 
    Joined -> True]}]

Here are some example outputs:

With NIID data Mathematica graphics

With AR(2) data Mathematica graphics

With a random walk data Mathematica graphics

share|improve this answer

Does the following direct implementation of your expressions work for you?

  ClearAll[data, dif, dbar, var, variogram];
  dif[list_] :=  Table[(list[[k + 1 ;;]] - list[[;; -(k+1)]]), 
  {k, Length[list] - 1}]
  var[list_] := Variance /@ (Most@(dif[list]));
  variogram[list_] := (#/First@#) &@var[list];

usage

 data = Accumulate[RandomReal[{-1, 1}, 100]];
 GraphicsRow[{ListLinePlot[data], ListLinePlot[variogram[data]]}]

enter image description here

Stationary data:

 data2 = RandomReal[{-1, 1}, 100];
 GraphicsRow[{ListLinePlot[data2], 
 ListLinePlot[Reverse@variogram[data2]]}]

enter image description here

share|improve this answer
    
This is not quite right. I don't think Differences works correctly for this application for differences bigger than 1. –  Mr Alpha Feb 28 '12 at 21:18
    
@MrAlpha, you are right. I changed the function to calculate the differences. Hope this one works. –  kguler Feb 28 '12 at 23:00

My literal implementation taken directly from the equations is below:

ClearAll[gk];
gk[z_,k_]:=Module[{n,s2k,dk,s21,d1},
  n = Length@z;
  dk[t_] := z[[t+k]]-z[[t]];
  s2k = Sum[(dk[t]-Sum[dk[t],{t,n-k}]/(n-k))^2,{t,n-k}]/(n-k-1);
  d1[t_] := z[[t+1]]-z[[t]];
  s21 = Sum[(d1[t]-Sum[d1[t],{t,n-1}]/(n-1))^2,{t,n-1}]/(n-2);
  s2k/s21
]

and this will generate an example plot

With[{z = RandomReal[{0, 1}, 100]},
  ListLinePlot[{gk[z, #] & /@ Range[98], z}]
]

enter image description here

Note that in the third equation, I am assuming the sum is from 1 to n - k.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.