/--------------------------------------------------------------------------------
This code is to implement MCMC inference on the models introduced in the paper

"A New Class of Nonseparable Space-time Covariance Models"

by Thais Fonseca and Mark Steel, 
University of Warwick, CRiSM Working Paper 08-13, 2008.

The wind-speed data as used in the paper are also included.
/--------------------------------------------------------------------------------

/--------------------------------------------------------------------------------
This code is written in OX that is an object-oriented statistical system. Ox comes 
in two versions: Ox Professional and Ox Console. Ox is available for Windows, Linux, 
Mac (OS X), and several Unix platforms. The latest version is 5.00 and the homepage 
is www.doornik.com. In the webpage you can find online documentation and the download 
links. The download is free for academics and available for a 2 week trial period to non-academics. 

The code AsymCovMod1.ox runs the model with generalized matern covariance in space.

The code AsymCovMod2and3.ox runs the matern model (model==2) or the cauchy model (model==3).

Both codes run the symmetric model (eps=0) or the asymmetric model (eps>0). And they 
also run the model with zero mean (EstimateDelta==0) or the model with covariates (EstimateDelta==1).
/--------------------------------------------------------------------------------
 

/--------------------------------------------------------------------------------
/******************** Functions used in the code *******************************/
/--------------------------------------------------------------------------------

/----------------------------------------------------------------/
/******************* Calculate distance in space ****************/
/----------------------------------------------------------------/
DistS(s,I)

This function calculates the Norm of s (matrix with I rows and 2 columns). I is the 
number of locations in space. 

/----------------------------------------------------------------/
/******** Calculate distance in space Gamma_1(s-epst) ***********/
/----------------------------------------------------------------/
DistSasym(s,I,t,eps)

This function calculates the Norm of (s_1-eps t,s_2), where eps is a scalar parameter
 and t is the time point.

/----------------------------------------------------------------/
/******************* Calculate distance in time *****************/
/----------------------------------------------------------------/
DistT(t,J)

This function calculates the distances in time. t is a vector with J rows.

/----------------------------------------------------------------/
/****************** Calculate Cauchy Covariance Function ********/
/----------------------------------------------------------------/
M1cauchy(space,lambda1,a,alpha)

This function calculates the Cauchy covariance function: (1 + ((space/a)^alpha)/a1) .^ (-lambda1).

/----------------------------------------------------------------/
/****************** Calculate Matern Covariance Function ********/
/----------------------------------------------------------------/
M1matern(space,lambda1,a,alpha)

This function calculates the Matern covariance function: 1/(2^(lambda1-1)*gammafact(lambda1 * 
(4*a1*(space/a) .^ alpha) .^ (lambda1/2) .* bessel((4*a1*(space/a) .^ alpha) .^ 0.5,lambda1) 

/----------------------------------------------------------------/
/****** Calculate Generalized Matern Covariance Function ********/
/----------------------------------------------------------------/
M1gig(space,lambda1,a,w,alpha

This function calculates the Generalized Matern covariance function: (1+1/w*(space/a) .^ alpha) 
.^ (-lambda1/2) .* bessel((4*w*(w+(space/a) .^ alpha)) .^ 0.5,lambda1)/bessel(2*w,lambda1) ;

/----------------------------------------------------------------/
/******************* Calculate the Likelihood Function **********/
/----------------------------------------------------------------/
likelihood(sig2,a,alpha,b,beta,vecz,Lag,tau2,lamb0,lamb1,w,eps,locat)

locat are the locations in space: matrix I by 2


/--------------------------------------------------------------------------/
/************************ VERY IMPORTANT **********************************/
/--------------------------------------------------------------------------/

/----------------------------------------------------------------/
/*************************** INPUTS *****************************/
/----------------------------------------------------------------/

  //// reading data
  vecz = loadmat("vecz.mat",1);     ////// VECTOR WITH COMPONENTS VEC(Z)
  locat = loadmat("s.mat");         ////// LOCATIONS IN SPACE I BY 2	  
  // times = loadmat("t.mat",1);    ////// IF IT IS THE CASE
  X = loadmat("X.mat");             ////// MATRIX WITH COVARIATES (INCLUDING FIRST COLUM OF ONES)
  Dhat = loadmat("DeltaHat.mat",1); ////// AN ESTIMATE FOR THE COEFICIENTS IN THE MEAN STRUCTURE (CAN BE ALL ZEROS)

// The files s.mat, X.mat and DeltaHat.mat need to have the number of rows and columns in the first row of the file 

  //// inputs
  I = 11; J = 10*365; N = I*J; M = 100000; 
  mylag = 3; EstimateDelta = 0; // if EstimaDelta=0 then the mean is zero
  model = 2; /// if model 2 then Matern if model 3 then cauchy

  /// initializing the parameters
  wk = 0.01;
  alphak = 1.6301;
  lamb2k = 1;  
  sigk = 0.35*ones(I,1);
  lamb1k = -1; thmcmc[0][1] = lamb1k;	
  tau2k = 0.10328; thmcmc[0][2] = tau2k;
  ak = 57.503; thmcmc[0][3] = ak;
  betak = 0.61619; thmcmc[0][4] = betak;
  bk = 2.1306; thmcmc[0][5] = bk;
  lamb0k = 1.1365; thmcmc[0][6] = lamb0k;
  if (EstimateDelta==0){  muk = zeros(N,1);} else{ deltak = Dhat[range(1,K,1)]; muk = X*deltak;}

  epsk = 1; eps[0] = epsk;

 
 
/----------------------------------------------------------------/
/*************************** OUTPUTS *****************************/
/----------------------------------------------------------------/
  
  savemat("Eps.mat",eps,1);                              ////// VECTOR WITH CHAIN OF EPS (NOW ZETA, M BY 1)
  savemat("DeltaMCMC.mat",delta,1);                      ////// MATRIX WITH CHAIN OF DELTA (M BY K)
  savemat("Sig2.mat",sig2,1);                            ////// MATRIX WITH CHAIN OF SIGMA2 (M BY I)
  savemat("th.mat",thmcmc,1);                            ////// MATRIX WITH CHAIN OF THETA (M BY 9)   												   //THETA=ALPHA,LAMBDA1,TAU2,A,BETA,B,LAMBDA0,LOGLIKELIHOOD,W)
  savemat("Accep.mat",(cont1/M)|(cont2/M)|(cont0/M),1);  ////// VECTOR WITH ACCEPTANCE RATE FOR EACH BLOCK