R/logpdfOU.R

In timedeppar: Infer Constant and Stochastic, Time-Dependent Model Parameters

#' Calculate log pdf of an Ornstein-Uhlenbeck process
#' 
#' This function calculates the log pdf of a realization of an Ornstein-Uhlenbeck process
#' with given parameters. 
#' The calculation can be done for an Ornstein-Uhlenbeck process with a random start value,
#' for a process conditional on the start value, or for a process conditional on start and end values.
#' The function includes the option of performing the calculation for lognormal marginal generated by 
#' exponential tranformation.
#' 
#' @param  t      vector of time points at which the OU process is available.
#' @param  y      vector of y-values corresponing to the t-values
#'                (note that t and y need to be of the same length).
#' @param  mean   asymptotic mean of the process.
#' @param  sd     asymptotic standard deviation of the process.
#' @param  gamma  rate coefficient for return to the mean
#' @param  cond   conditioning: 0 indicates no conditioning, 1 conditioning to start value, and
#'                2 conditioning to start and end value
#' @param  log    if true, the log pdf of the log of y is calculated
#'                (mean and sd are interpreted in y, not in log(y) units)
#' 
#' @return the function returns the log pdf
#' 
#' @examples 
#' OU <- randOU(mean=0,sd=1,gamma=1,t=0:1000/1000)
#' logpdfOU(OU$t,OU$y,mean=0,sd=1,gamma=1)


logpdfOU <- function(t,y,mean=0,sd=1,gamma=1,cond=0,log=FALSE)
{
  # consistency checks:
  
  n <- length(t)
  if ( length(y) != n )
  {
    warning("length(y) must be equal to length(t)")
    return(NA)
  }
  if ( min(diff(t)) <= 0 )
  {
    warning("t must be strictly increasing")
    return(NA)
  }
  if ( sd <= 0 )
  {
    warning("sd must be positive")
    return(NA)
  }
  if ( gamma <= 0 )
  {
    warning("gamma must be positive")
    return(NA)
  }
  
  logpdf <- 0
  if ( log ) 
  {
    # if log=TRUE calculate meanlog and sdlog (mean and sd are in original units),
    # call function with these parameters, and convert pdf to original units:
    
    if ( mean <= 0 )
    {
      warning("if log=TRUE, mean must be positive")
      return(NA)
    }
    if ( min(y) <= 0 )
    {
      warning("if log=TRUE, all y values must be positive")
      return(NA)
    }
    meanlog <- log(mean/(sqrt(1+sd*sd/(mean*mean))))
    sdlog   <- sqrt(log(1+sd*sd/(mean*mean)))
    logpdf  <- logpdfOU(t=t,y=log(y),mean=meanlog,sd=sdlog,gamma=gamma,
                        cond=cond,log=FALSE)
    logpdf <- logpdf - sum(log(y[2:(length(y)-1)]))
    if ( cond < 2 ) logpdf <- logpdf - log(y[length(y)])
    if ( cond < 1 ) logpdf <- logpdf - log(y[1])
  }
  else
  {  
    # if log=FALSE calculate log pdf for non-logarithmic case:
    
    if ( cond == 0 ) logpdf <- dnorm(y[1],mean=mean,sd=sd,log=TRUE)  # initial point
    if ( n > 1 )
    {
      if ( cond < 2 )  # iterative conditioning on previous point
      {
        # didactical version:
        
        # for ( i in 2:n )
        # {
        #   logpdf <- logpdf + dnorm(x    = y[i],
        #                            mean = mean+(y[i-1]-mean)*exp(-gamma*(t[i]-t[i-1])),
        #                            sd   = sd*sqrt(1-exp(-2*gamma*(t[i]-t[i-1]))),
        #                            log  = TRUE)
        # }
        
        # more efficient version:
        
        logpdf <- logpdf + sum(dnorm(x    = y[-1],
                                     mean = mean+(y[-length(y)]-mean)*exp(-gamma*diff(t)),
                                     sd   = sd*sqrt(1-exp(-2*gamma*diff(t))),
                                     log  = TRUE))
      }
      else
      {
        if ( n > 2 )  # iterative conditioning on previous point and end point
        {
          # didactical version:
          
          # for ( i in 2:(n-1) )
          # {
          #   logpdf <- logpdf + 
          #             dnorm(x    = y[i],
          #                   mean = mean+(y[i-1]-mean)*exp(-gamma*(t[i]-t[i-1]))*(1-exp(-2*gamma*(t[n]-t[i])))/
          #                                            (1-exp(-2*gamma*(t[n]-t[i-1]))) +
          #                               (y[n]-mean)  *exp(-gamma*(t[n]-t[i]))*(1-exp(-2*gamma*(t[i]-t[i-1])))/
          #                                             (1-exp(-2*gamma*(t[n]-t[i-1]))),
          #                   sd   = sd*sqrt((1-exp(-2*gamma*(t[i]-t[i-1])))*(1-exp(-2*gamma*(t[n]-t[i])))/
          #                                  (1-exp(-2*gamma*(t[n]-t[i-1])))),
          #                   log  = TRUE)
          # }
          
          # more efficient version:
          
          logpdf <- logpdf +
                    sum(dnorm(x    = y[-1],
                              mean = mean+(y[-length(y)]-mean)*exp(-gamma*diff(t))*(1-exp(-2*gamma*(t[n]-t[-1])))/
                                                               (1-exp(-2*gamma*(t[n]-t[-length(t)]))) +
                                          (y[n]-mean)         *exp(-gamma*(t[n]-t[-1]))*(1-exp(-2*gamma*diff(t)))/
                                                               (1-exp(-2*gamma*(t[n]-t[-length(t)]))),
                              sd   = sd*sqrt((1-exp(-2*gamma*diff(t)))*(1-exp(-2*gamma*(t[n]-t[-1])))/
                                             (1-2*gamma*(t[n]-t[-length(t)]))),
                              log  = TRUE))
            
        }
      }
    }
  }
  return(logpdf)
}