TraME-R: Transportation Methods for Econometrics

################################################################################
##
##   Copyright (C) 2015 - 2016 Alfred Galichon
##
##   This file is part of the R package TraME.
##
##   The R package TraME is free software: you can redistribute it and/or modify
##   it under the terms of the GNU General Public License as published by
##   the Free Software Foundation, either version 2 of the License, or
##   (at your option) any later version.
##
##   The R package TraME is distributed in the hope that it will be useful,
##   but WITHOUT ANY WARRANTY; without even the implied warranty of
##   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##   GNU General Public License for more details.
##
##   You should have received a copy of the GNU General Public License
##   along with TraME. If not, see <http://www.gnu.org/licenses/>.
##
################################################################################
#
# References:
#
# A. Galichon, S. Weber: "Estimation of Matching Function Equilibria"
#           
#           
# mmfs has n, and m and neededNorm
# M.mmfs <- function(mmfs, ax, by)
#
################################################################################
########################    Default and generic methods ########################
################################################################################
#
margxInv.default <- function(xs, mmfs, theMu0ys)
{
  nbX = length(mmfs$n)
  if (is.null(mmfs$neededNorm))
  {
    coeff = 1
    ubs = mmfs$n
  }
  else
  {
    coeff = 0
    ubs = rep(1e10,nbX)
  }
  #
  if(is.null(xs)){
    xs = 1:nbX
  }
  theaxs = rep(0,length(xs))
  #
  i = 0
  for(x in xs){
    i = i+1
    root_fn <- function(z) (ifelse(coeff,z,0) - mmfs$n[x] + sum(M(mmfs,z,theMu0ys,xs=x)))
    theaxs[i] = uniroot(root_fn, c(0,ubs[x]), tol = 1e-300)$root # Keith: fix tolerence   
  }
  #
  return(theaxs)
}
#
margyInv.default <- function(ys, mmfs, theMux0s)
{
  nbY = length(mmfs$m)
  if (is.null(mmfs$neededNorm))
  {
    coeff = 1
    ubs = mmfs$m
  }
  else
  {
    coeff = 0
    ubs = rep(1e10,nbY)
  }
  #
  if(is.null(ys)){
    ys = 1:nbY
  }
  thebys = rep(0,length(ys))
  #
  j = 0
  for(y in ys){
    j = j+1
    root_fn <- function(z) (ifelse(coeff,z,0) - mmfs$m[y] + sum(M(mmfs,theMux0s,z,ys=y)))
    thebys[j] = uniroot(root_fn, c(0,ubs[y]), tol=1e-300)$root
  }
  #
  return(thebys)
}
#
################################################################################
########################           geo MMfs             ########################
########################  (generated by TU transfers)   ########################
################################################################################
build_geommfs <- function(n,m,K,neededNorm=NULL)
{
  ret = list(n=n,
             m=m,
             nbX=length(n),
             nbY=length(m),
             neededNorm=neededNorm,
             K = K
  )
  class(ret)="geommfs"
  return(ret)
}
#
mmfsTranspose.geommfs <- function(mmfs)
{
  ret = list(n=mmfs$m,
             m=mmfs$n,
             nbX = mmfs$nbY,
             nbY = mmfs$nbX,
             neededNorm=normalizationTranspose(mmfs$neededNorm),
             K = t(mmfs$K)
  )
  class(ret)="geommfs"
  return(ret)
}
#
M.geommfs <- function(mmfs, mux0s, mu0ys, xs=1:length(mmfs$n), ys=1:length(mmfs$m))
{
  term_1 = mmfs$K[xs,ys]
  term_2 = sqrt(mux0s %*% t(mu0ys))
  #
  ret = term_1 * term_2
  #
  return(ret)
}
# 
dmux0s_M.geommfs <- function (mmfs,mux0s,mu0ys)
{
  term_1 = mmfs$K / 2
  term_2 = sqrt( (1/mux0s) %*% t(mu0ys))
  #
  ret = term_1 * term_2
  #
  return(ret)
}
#
dmu0ys_M.geommfs <- function (mmfs,mux0s,mu0ys)
{
  term_1 = mmfs$K / 2
  term_2 = sqrt( mux0s %*% t(1/mu0ys))
  #
  ret = term_1 * term_2
  #
  return(ret)
}
#
dparams_M.geommfs <- function (mmfs,mux0s,mu0ys,deltaparamsM=NULL)
{
  ret <- 0
  if(is.null(deltaparamsM)){
    ret = Diagonal(sqrt(mux0s %*% t(mu0ys)),n=mmfs$nbX*mmfs$nbY)
  }else{
    ret = c(matrix(deltaparamsM,nrow = nbX) * sqrt(mux0s %*% t(mu0ys)) )
  }
  return(ret)
}
#
margxInv.geommfs <- function(xs, mmfs, theMu0ys)
{
  if (is.null(xs)) {xs = 1:length(mmfs$n)}
  #
  sqrtBys = sqrt(theMu0ys)
  if(is.null(mmfs$neededNorm)){
    b = (mmfs$K[xs,] %*% sqrtBys)/2
    sqrtAxs = sqrt(mmfs$n[xs]+ b*b) - b
  }else{
    sqrtAxs = mmfs$n / c(mmfs$K[xs,] %*% sqrtBys)
  }
  #
  ret = c(sqrtAxs*sqrtAxs)
  #
  return(ret)
}
#
margyInv.geommfs <- function(ys, mmfs, theMux0s)
{
  if (is.null(ys)) {xs = 1:length(mmfs$m)}
  #
  sqrtAxs = sqrt(theMux0s)
  if(is.null(mmfs$neededNorm)){
    b = c(sqrtAxs %*% mmfs$K[,ys])/2
    sqrtBys = sqrt(mmfs$m[ys] + b*b) - b
  }else{
    sqrtBys = mmfs$m / c(sqrtAxs %*% mmfs$K[,ys])
  }
  #
  ret = sqrtBys*sqrtBys
  #
  return(ret)
}

################################################################################
########################            min MMfs            ########################
########################  (generated by NTU transfers)  ########################
################################################################################
build_minmmfs <- function(n,m,A,B,neededNorm=NULL)
{
  ret = list(n=n,
             m=m,
             nbX=length(n),
             nbY=length(m),
             neededNorm=neededNorm,
             A = A,
             B = B)
  class(ret)="minmmfs"
  return(ret)
}
#
mmfsTranspose.minmmfs <- function(mmfs)
{
  ret = list(n = mmfs$m,
             m = mmfs$n,
             nbX = mmfs$nbY,
             nbY = mmfs$nbX,
             neededNorm = normalizationTranspose(mmfs$neededNorm),
             A = t(mmfs$B),
             B = t(mmfs$A)
  )
  class(ret)="minmmfs"
  return(ret)
}
#
M.minmmfs <- function(mmfs, mux0s, mu0ys, xs=1:length(mmfs$n), ys=1:length(mmfs$m))
{
  term_1 = mux0s * mmfs$A[xs,ys]
  term_2 = t( mu0ys * t(mmfs$B[xs,ys] ))
  #
  ret = pmin(term_1, term_2)
  #
  return(ret)
}
# 
dmux0s_M.minmmfs <- function (mmfs,mux0s,mu0ys)
{
  term_1 = mux0s * mmfs$A
  term_2 = t( mu0ys * t(mmfs$B ))
  #
  return(ifelse(term_1 <= term_2,1,0) * mmfs$A)
}
#
dmu0ys_M.minmmfs <- function (mmfs,mux0s,mu0ys)
{
  term_1 = mux0s * mmfs$A
  term_2 = t( mu0ys * t(mmfs$B ))
  #
  return(ifelse(term_1 >= term_2,1,0) * mmfs$B)  
}
#
dparams_M.minmmfs <- function (mmfs,mux0s,mu0ys,deltaparamsM=NULL)
{
  term_1 = mux0s * mmfs$A
  term_2 = t( mu0ys * t(mmfs$B ))
  t1lessthant2 = ifelse(term_1 <= term_2,1,0)
  der1 = mux0s * t1lessthant2
  der2 = t(mu0ys * t(1-t1lessthant2) )
  ret <- 0
  if(is.null(deltaparamsM)){
    ret = cbind(Diagonal(x=der1),
                Diagonal(x=der2 ))
  }else{
    deltaparams1 = matrix(deltaparamsM[1:(mmfs$nbX*mmfs$nbY)],nrow = nbX)
    deltaparams2 = matrix(deltaparamsM[(1+mmfs$nbX*mmfs$nbY):(2*mmfs$nbX*mmfs$nbY)],nrow = nbX)
    ret = c(deltaparams1 * der1 + deltaparams2 * der2 )
  }
  return(ret)
}
#
margxInv.minmmfs = function(xs,mmfs,theMu0ys)
{
  if (is.null(xs)) {xs = 1:length(mmfs$n)}
  if (!is.null(mmfs$neededNorm)) {stop('not supported yet')}
  theaxs = inversePWA(mmfs$n[xs], t ( t(mmfs$B[xs,] / mmfs$A[xs,])* theMu0ys ),mmfs$A[xs,]) 
  return(theaxs)
}
#
margyInv.minmmfs = function(ys,mmfs,theMux0s)
{
  if (is.null(ys)) {ys = 1:length(mmfs$m)}
  if (!is.null(mmfs$neededNorm)) {stop('not supported yet')}
  thebys = inversePWA(mmfs$m[ys], t(( mmfs$A[,ys] / mmfs$B[,ys]) * theMux0s), t(mmfs$B[,ys] ) ) 
  return(thebys)
}

################################################################################
########################            cod MMfs            ########################
########################  (generated by LTU transfers)  ########################
################################################################################
build_codmmfs <- function(n,m,lambda,K,neededNorm=NULL)
{
  ret = list(n=n,
             m=m,
             nbX=length(n),
             nbY=length(m),
             neededNorm=neededNorm,
             lambda = lambda,
             K=K,
             aux_zeta = 1-lambda)
  class(ret)="codmmfs"
  return(ret)
}
#
mmfsTranspose.codmmfs <- function(mmfs)
{
  ret = list(n=mmfs$m,
             m=mmfs$n,
             nbX = mmfs$nbY,
             nbY = mmfs$nbX,
             neededNorm=normalizationTranspose(mmfs$neededNorm),
             lambda = t(mmfs$aux_zeta),
             K=t(mmfs$K),
             aux_zeta = t(mmfs$lambda)
  )
  class(ret)="codmmfs"
  return(ret)
}
#
M.codmmfs <- function(mmfs, mux0s, mu0ys, xs=1:length(mmfs$n), ys=1:length(mmfs$m))
{
  term_1 = mux0s^mmfs$lambda[xs,ys]
  term_2 = t( mu0ys^t(mmfs$aux_zeta[xs,ys]) )
  term_3 = mmfs$K[xs,ys]
  #
  ret = term_1 * term_2 * term_3
  #
  return(ret)
}
#
# 
dmux0s_M.codmmfs <- function (mmfs,mux0s,mu0ys)
{
  term_1 = mmfs$lambda* mux0s^(mmfs$lambda-1)
  term_2 = t( mu0ys^t(mmfs$aux_zeta) )
  term_3 = mmfs$K
  #
  ret = term_1 * term_2 * term_3
  #
  return(ret)
}
#
dmu0ys_M.codmmfs <- function (mmfs,mux0s,mu0ys)
{
  term_1 = mux0s^mmfs$lambda
  term_2 = mmfs$aux_zeta * t( mu0ys^t( mmfs$aux_zeta -1) )
  term_3 = mmfs$K
  #
  ret = term_1 * term_2 * term_3
  #
  return(ret)
}
#
dparams_M.codmmfs <- function (mmfs,mux0s,mu0ys,deltaparamsM=NULL)
{
  term_1 = mux0s^(mmfs$lambda)
  term_2 = t( mu0ys^t(mmfs$aux_zeta) )
  term_3 = mmfs$K
  logratio = log( mux0s %*% t(1 / mu0ys))
  der1 = logratio * term_1*term_2*term_3
  der2 = term_1 * term_2
  
  ret <- 0
  if(is.null(deltaparamsM)){
    ret = cbind (Diagonal(x = der1),
                 Diagonal(x = der2))
  }else{
    deltaparams1 = matrix(deltaparamsM[1:(mmfs$nbX*mmfs$nbY)],nrow = nbX)
    deltaparams2 = matrix(deltaparamsM[(1+mmfs$nbX*mmfs$nbY):(2*mmfs$nbX*mmfs$nbY)],nrow = nbX)
    ret = c(deltaparams1 * der1 + deltaparams2 * der2)}
    
    return(ret)
}

#
################################################################################
########################            ces MMfs            ########################
################################################################################
# corresponds to ETU transfers
build_cesmmfs <- function(n,m,C,D,kappa,neededNorm=NULL)
{
  ret = list(n=n,
             m=m,
             nbX=length(n),
             nbY=length(m),
             neededNorm=neededNorm,
             C = C,
             D = D,
             kappa = kappa # note that kappa is -1/tau in the ETU-logit interpretation
  )
  class(ret)="cesmmfs"
  return(ret)
}
#
mmfsTranspose.cesmmfs <- function(mmfs)
{
  ret = list(n = mmfs$m,
             m = mmfs$n,
             nbX = mmfs$nbY,
             nbY = mmfs$nbX,
             neededNorm = normalizationTranspose(mmfs$neededNorm),
             C = t(mmfs$D),
             D =  t(mmfs$C),
             tauinv = t(mmfs$tauinv)
  )
  class(ret)="cesmmfs"
  return(ret)
}
#
M.cesmmfs <- function(mmfs, mux0s, mu0ys, xs=1:length(mmfs$n), ys=1:length(mmfs$m))
{
  term_1 = mmfs$C[xs,ys] * (mux0s^mmfs$kappa[xs,ys])
  term_2 = mmfs$D[xs,ys] * t(mu0ys^t(mmfs$kappa[xs,ys]))
  #
  ret = ((term_1 + term_2)/2)^(1/mmfs$kappa[xs,ys])
  #
  return(ret)
}
#
dmux0s_M.cesmmfs <- function (mmfs,mux0s,mu0ys)
{
  term_1 = mmfs$C * (mux0s^mmfs$kappa)
  term_2 = mmfs$D * t(mu0ys^t(mmfs$kappa))
  #
  mus = ((term_1 + term_2)/2)^(1/mmfs$kappa)
  #
  
  return( mmfs$C * (mus / mux0s)^(1 - mmfs$kappa) /2 )
}
#
dmu0ys_M.cesmmfs <- function (mmfs,mux0s,mu0ys)
{
  term_1 = mmfs$C * (mux0s^mmfs$kappa)
  term_2 = mmfs$D * t(mu0ys^t(mmfs$kappa))
  #
  mus = ((term_1 + term_2)/2)^(1/mmfs$kappa)
  #
  
  return( mmfs$D * t( (t(mus) / mu0ys)^t(1 - mmfs$kappa) ) /2 )
}
#

dparams_M.cesmmfs <- function (mmfs,mux0s,mu0ys,deltaparamsM=NULL)
{
  term_1 = mmfs$C * (mux0s^mmfs$kappa)
  term_2 = mmfs$D * t(mu0ys^t(mmfs$kappa))
  #
  mus = ((term_1 + term_2)/2)^(1/mmfs$kappa)
  
  num = (mus / mmfs$kappa) * ( mmfs$C * mux0s^mmfs$kappa * log(mux0s) + mmfs$D * t(mu0ys^t(mmfs$kappa)*log(mu0ys)) - log(mus) )
  denom = mmfs$C * mux0s^mmfs$kappa  + mmfs$D * t(mu0ys^t(mmfs$kappa))
  
  der1 = mus^(1 - mmfs$kappa) * mux0s^mmfs$kappa /2
  der2 = mus^(1 - mmfs$kappa) * t(mu0ys^t(mmfs$kappa))/ 2
  der3 = num / denom
  
  ret <- 0
  if(is.null(deltaparamsM)){
    ret = cbind (Diagonal(x = der1),
                 Diagonal(x = der2),
                 Diagonal(x = der3 ))
  }else{
    deltaparams1 = matrix(deltaparamsM[1:(mmfs$nbX*mmfs$nbY)],nrow = nbX)
    deltaparams2 = matrix(deltaparamsM[(1+mmfs$nbX*mmfs$nbY):(2*mmfs$nbX*mmfs$nbY)],nrow = nbX)
    deltaparams3 = matrix(deltaparamsM[(1+2*mmfs$nbX*mmfs$nbY):(3*mmfs$nbX*mmfs$nbY)],nrow = nbX)
    
    ret = c(deltaparams1 * der1 + deltaparams2 * der2 + deltaparams3 * der3 )
  }
  
  return(ret)
}