// Flaremodel.
#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() ()
{
  
  using namespace density;
  
  // Data input
  DATA_MATRIX(y);                                                     // Flaring data
  DATA_IMATRIX(censor);                                               // Censor data with 1 for non-censored, 0 for cencored
  DATA_VECTOR(price);                                                 // Price covariate
  
  // Parameters
  PARAMETER(alpha1);                                                  // effect of price
  PARAMETER(log_sigma);       Type sigma = exp(log_sigma);            // sd of residual 
  PARAMETER(log_sigma_alpha); Type sigma_alpha = exp(log_sigma_alpha);// sd of random effect (alpha0)
  PARAMETER_VECTOR(alpha0);                                           // random effect
  
  // Get SD's of:
  ADREPORT(alpha1);
  ADREPORT(sigma);
  ADREPORT(sigma_alpha);
  
  
  int n = y.rows();        // number of fields
  int time = y.cols();  // number of years

  
  Type nll = 0.0;          // Initiate negative log-likelihood
  
  for(int i=0; i< n; i++){
    for(int t = 0; t < time; t++){
      
      if (y(i,t) != y(i,t)) { // f != f returns true if and only if f is NaN. 
        
        // Replace missing values (NA in R, NaN in C++) with 0 contribution in likelihood
        nll -= 0;
      }
      else { 
        Type eta = alpha0(i) + alpha1*price(t);
        nll -= censor(i,t)*dnorm(Type(y(i,t)), eta, sigma, true) + (1 - censor(i,t))*log(pnorm(Type(0), eta, sigma));
            } // end of if-else
    } // end of looping over t
    nll -= dnorm(Type(alpha0(i)), Type(0), sigma_alpha, true);
  } // end of loop over i
  
 REPORT(alpha0);
  
  return nll;
}