content: module QuantLib.Math.Copulas
          ( Copula (..),
            Copulas (..)
          ) where
  
  {-| Copula type class.
   -| Normally instance should implement only copulaFunc.
   -| Method copula provides a precheck for [0..1] range for x and y but real implementation is in copulaFunc
   -}
  class Copula a where
          copula :: a -> Double -> Double -> Maybe Double
          copula t = precheckRange (copulaFunc t)
          copulaFunc :: a -> Double -> Double -> Maybe Double
  
  -- Copula must be in [0,1] range
  precheckRange :: (Double->Double->Maybe Double) -> Double -> Double -> Maybe Double
  precheckRange f x y
          | x >= 0.0 || x <= 1.0
          || y>= 0.0 || y <= 1.0
                  = f x y
          | otherwise
                  = Nothing
  
  {-| Copula data types with parameters required by the concrete copula definition
   -}
  data Copulas = ClaytonCopula Double -- ^ Clayton copula
          | MinCopula -- ^ Min copula
          | MaxCopula -- ^ Max copula
          | AliMikhailHaqCopula Double -- ^ Ali-Mikhail-Haq copula
          | FarlieGumbelMorgensternCopula Double -- ^ Farlie-Gumbel-Morgenstern copula
          | FrankCopula Double -- ^ Frank copula
          | GalambosCopula Double -- ^ Galambos copula
          | GaussianCopula Double -- ^ Gaussian copula /Not implemented yet!/
          | GumbelCopula Double -- ^ Gumbel copula
          | HuslerReissCopula Double -- ^ Husler-Reiss copula /Not implemented yet!/
          | IndependentCopula -- ^ Independent copula
          | MarshallOlkinCopula Double Double -- ^ Marshall-Olkin copula
          | PlackettCopula Double -- ^ Plackett copula
  
  instance Copula Copulas where
          copulaFunc (ClaytonCopula theta)                    = claytonCopula theta
          copulaFunc MinCopula                                = minCopula
          copulaFunc MaxCopula                                = maxCopula
          copulaFunc (AliMikhailHaqCopula theta)              = aliMikhailHaqCopula theta
          copulaFunc (FarlieGumbelMorgensternCopula theta)    = farlieGumbelMorgenstern theta
          copulaFunc (FrankCopula theta)                      = frankCopula theta
          copulaFunc (GalambosCopula theta)                   = galambosCopula theta
          copulaFunc (GaussianCopula rho)                     = gaussianCopula rho
          copulaFunc (GumbelCopula theta)                     = gumbelCopula theta
          copulaFunc (HuslerReissCopula theta)                = huslerReissCopula theta
          copulaFunc IndependentCopula                        = independentCopula
          copulaFunc (MarshallOlkinCopula a b)                = marshallOlkinCopula a b
          copulaFunc (PlackettCopula theta)                   = plackettCopula theta
  
  {- Private implementations   -}
  
  {-# ANN aliMikhailHaqCopula "NoHerbie" #-}
  aliMikhailHaqCopula :: (Fractional a, Ord a) => a -> a -> a -> Maybe a
  aliMikhailHaqCopula theta x y
          | theta >= -1.0 && theta <= 1.0
                  = Just ((x*y)/(1.0 - theta * (1.0-x)*(1.0-y)))
          | otherwise
                  = Nothing
  
  {-# ANN farlieGumbelMorgenstern "NoHerbie" #-}
  farlieGumbelMorgenstern :: (Fractional a, Ord a) => a -> a -> a -> Maybe a
  farlieGumbelMorgenstern theta x y
          | theta >= -1.0 && theta <= 1.0
                  = Just (x*y + theta*x*y*(1.0-x)*(1.0-y))
          | otherwise
                  = Nothing
  
  {-|  Original code and algorithm from the Quantlib project
       implemented in Haskell by Nicholas Pezolano
                                    npezolano "at" gmail.com
  -}
  {-# ANN claytonCopula "NoHerbie" #-}
  claytonCopula :: Double -> Double -> Double -> Maybe Double
  claytonCopula theta x y
    |  theta ==0
    || theta < -1.0
     = Nothing
  
     | otherwise
     = Just $ max( (x ** (-theta) ) +   (y ** (-theta)-1.0)  **   (-1.0/theta)) 0
  
  minCopula ::  Ord a => a -> a -> Maybe a
  minCopula x y = Just (min x y)
  
  maxCopula ::  (Fractional a, Ord a) => a -> a -> Maybe a
  maxCopula x y = Just (max 0.0 (x+y-1.0))
  
  frankCopula ::  (Eq a, Floating a) => a -> a -> a -> Maybe a
  frankCopula theta x y
      | theta     == 0.0 = Nothing
      | otherwise = Just (-1.0/theta * log (1 + (exp (-theta*x) - 1.0) * (exp (-theta*y) -1.0) / (exp (-theta) - 1.0)   ))
  
  {-# ANN galambosCopula "NoHerbie" #-}
  galambosCopula ::  (Floating a, Ord a) => a -> a -> a -> Maybe a
  galambosCopula theta x y
      | theta <= 0.0  = Nothing
      | otherwise     = Just (x * y * exp ( ( (-log x) ** (-theta) + (-log y) ** (-theta) ) ** (-1.0/theta)  ))
  
  gaussianCopula ::  (Fractional a, Ord a) => a -> t -> t1 -> Maybe a1
  gaussianCopula rho _ _
      | rho >= -1.0 && rho <= 1.0 = undefined
      | otherwise                 = Nothing
  
  {-# ANN gumbelCopula "NoHerbie" #-}
  gumbelCopula ::  (Floating a, Ord a) => a -> a -> a -> Maybe a
  gumbelCopula theta x y
      | theta >= 1.0  = Just (exp ( - ( (-log x) ** theta + (-log y) ** theta) ** (1.0/theta)))
      | otherwise     = Nothing
  
  huslerReissCopula :: (Fractional a, Ord a) => a -> a -> a -> Maybe a
  huslerReissCopula theta _ _
      | theta > 0.0   = undefined
      | otherwise     = Nothing
  
  independentCopula ::  Num a => a -> a -> Maybe a
  independentCopula x y = Just (x*y)
  
  marshallOlkinCopula :: (Floating a, Ord a) => a -> a -> a -> a -> Maybe a
  marshallOlkinCopula a b x y
      | a >= 0.0 && b >= 0.0  = Just (min (y * (x ** (1-a))) (x * (y ** (1-b))))
      | otherwise             = Nothing
  
  {-# ANN plackettCopula "NoHerbie" #-}
  plackettCopula ::  (Floating a, Ord a) => a -> a -> a -> Maybe a
  plackettCopula theta x y
      | theta >= 0.0 && theta /= 1.0  = Just $ (sumXyTheta1 - sqrt (sumXyTheta1 * sumXyTheta1 - 4.0 * x * y * theta * theta1))/(2*theta1)
      | otherwise                     = Nothing
          where   sumXy           = x + y
                  theta1          = theta - 1.0
                  sumXyTheta1     = 1.0 + theta1 * sumXy
  