/*============================================================================

 Function      ar1

 Usage         ar1(nz, lambda, mu, sigma, phi, Z, P)

 Arguments     nz:     constant integer representing number of values in the
                       AR1 approximation grid.

	       lambda: constant double representing the number of standard
	               deviations from the mean to place the bounds of the AR1
		       grid.
	           
	       mu:     constant double representing the mean of the AR1
	               process.

	       sigma:  constant double representing the standard deviation of
	               the AR1 process.

	       rho:    constant double representing the persistence of the AR1
	               process.

	       Z:      pointer to array of doubles which stores the AR1 grid
	               grid values.

	       P:      pointer to array of doubles which stores the AR1
	               transition matrix values.
	              
 Description   This function is a CUDA kernel that computes the grid values
               and transition matrix for a finite approximation of an AR1
	       process according to the method of Tauchen (1986).

 Return value  void.

 =============================================================================

 Author:       Eric M. Aldrich

 Contact:      ealdrich@gmail.com

 Date:         20 April 2010

 ============================================================================*/

__global__ void ar1(const int nz, const double lambda, const double mu,
		    const double sigma, const double rho, double* Z, double* P)
{ 

  // thread
  const int i = blockIdx.x * blockDim.x + threadIdx.x;
  int j;

  // grid for TFP
  const double sigma_z = sigma/sqrt(1-pow(rho,2));
  const double mu_z = mu*(1/(1-rho));
  const double zmin = mu_z - lambda*sigma_z;
  const double zmax = mu_z + lambda*sigma_z;
  const double zstep = (zmax-zmin)/(nz-1);
  Z[i] = exp(zmin + zstep*i);

  __syncthreads();

  // transition matrix
  double normarg1, normarg2;
  P[i*nz] = ncdf((zmin - mu - rho*log(Z[i]))/sigma + 0.5*zstep/sigma);
  P[i*nz+nz-1] = 1-P[i*nz];
  for(j = 1 ; j < nz-1 ; ++j){
    normarg1 = (log(Z[j]) - mu - rho*log(Z[i]))/sigma + 0.5*zstep/sigma;
    normarg2 = (log(Z[j]) - mu - rho*log(Z[i]))/sigma - 0.5*zstep/sigma;
    P[i*nz+j] = ncdf(normarg1) - ncdf(normarg2);
    P[i*nz+nz-1] -= P[i*nz+j];
  }
}

