Description of R codes for the paper "A test for stationarity for spatio-temporal data" by S. Bandyopadhyay, C. Jentsch and Suhasini Subba Rao (2016).

The code requires the prior installation and loading of the R packages `MASS' and `fields' (in R, go to Packages -> Install package(s)..., follow the instructions, and then issue the command library(MASS), library(fields)).

The code is deliberately kept simple so that it is hopefully easy to understand and use. Please do let us know (email: sob210@lehigh.edu) if you make any modifications or improvements to it or if you have any questions. 

Also please refer all use of the code. 

The main function is 'analysis.r' which the user should use.

################################################################
THIS 'README.TXT' GIVES DETAILED DESCRIPTION OF DIFFERENT FILES IN THE FOLDER.
################################################################

1. analysis.r: This file has the functions which reports all 10 test statistics along with the associated P-values for a given data as described in the paper. 

We have sourced all other functions in the folder in this file.


2. simdata.Rdata: This is a simulated dataset generated using model NS3 for n = 100, T = 200, and lambda = 5. 

This .Rdata file contains a list sim.data where:
(a) site = sim.data[[1]] gives a nx2 matrix of locations
(b) data = sim.data[[2]] gives a nxT matrix of spatio-temporal data.

3. R-grd-fn.r: The function gives r-grid values for the choice of spatial component in \mathcal{P} and \mathcal{P}^{'}. For example, in the code

r.grd[[1]] = {(1, 0), (1, 1), (0, 1), (-1, 1)}
r.grd[[2]] = {(2, 0), (2, 1), (2, 2), (1, 2), (0, 2), (-1, 2),
		  (-2, 2), (-2 1)}


4. a-g-hat-fn.r: The function is used for the calculation of a_{g}(\omega_{k},r_{1},r_{2}) for all k and for a particular g() function g(\Omega) = \sum_{j=1}^{L}\exp{i v_{j}^{'}\Omega} as described in the paper. The v grid for the calculation of g() function has been calculated in the analysis.r file. 

Note that the data comes here as fft.dat - each row of this data set is one site over time, so each column is the fourier transform at frequency omega_{k}, i.e.,

#\begin{eqnarray*}
#a_{g}(\omega_{k_{2}};\rb_{1},r_{2}) =
#\frac{1}{n^{2}}\sum_{j_{1},j_{2}=1}^{n} \prod_{s=1}^{d}D_{n}\left(\frac{s_{d,j_{1}} -  s_{d,j_{2}} - v_{d}}{\lambda}\right)
#J_{\lambda,\bs_{j_1}}(\omega_{k_{2}})
#\overline{J_{\lambda,\bs_{j_2}}(\omega_{k_{2}+r_{2}})}e^{-i\bs_{j_{2}}\ob_{\rb_{1}}}.
#\end{eqnarray*}

5. all-test-statistics-fn.r: In this function we calculate the 10 test statistics as reported in the paper. 

In our simulation study (and also in the example), we caluclate all test statistics without the nugget term.

6. crit-val-fn.r: The following function generates samples from the asymptotic ‘finite sample’ approximations of the distribution of the test statistics under the null. The values are needed for obtaining the P-value of the tests.