R/Examplehighdimdiffindiff.R

In psdsam/HDdiffindiff: High Dimensional Doubly Robust Diff-in-Diff Estimator

#' An example to implement Dr-did estimator
#' @param p dimension of the linear specification
#' @param q parameter for the nonparametric specification
#' @param n0 number of sample size
#' @param rho.X parameter to construct the correlation matrix for the high dimensional covariates
#' @param method basis function in nonparametric specification:
#' \itemize{
#' \item (default) Pol: Polynomial basis with parameter q/2 as the degree of polynomial;
#' \item Tri: Trigonometric polynomials with parameter q/2 as the degree of polynomial;
#' }
#' @export Examplehighdimdiffindiff

Examplehighdimdiffindiff <- function(p=50, q=8, n0=500, method = "Pol", rho.X=0.5){

  omega0 = rep(0,p)
  for (i in 1:10){
    omega0[i] = 1/i
  }
  omega0 = matrix(omega0, p, 1)
  omega1 = rep(0,p)
  for (i in 1:min(p, 15)){
    omega1[i] = 2/i
  }
  omega1 = matrix(omega1, p, 1)
  theta0 = rep(0,p)
  for (i in 1:10){
    theta0[i] = 1/i
  }

  z = rnorm(n0,0, 1)

  alpha = 0
  sig   <-  0.5^(p-toeplitz(p:1))
  x = mvrnorm(n0, mu = rep(0, p), Sigma = sig)
  y0 = rnorm(n0,0,1)*(1/sqrt(2)*z +1/sqrt(2)*x[,1])
  prop = 1-1/(1+exp( x %*% theta0))
  treat = rbinom(n0,1,prop)

  z = rnorm(n0,0, 1)
  epsilon0 = rnorm(n0,0,1)
  epsilon1 = rnorm(n0,0,2)

  Phi0 = (x %*% omega0)
  Phi1 = (x %*% omega1) +  exp(z)

  y1 = y0 + alpha + (Phi1 + epsilon1)*treat + (Phi0 + epsilon0)*(1-treat)


  ff = highdimdiffindiff_crossfit(y0, y1, treat, x, z,k=3, method =method, q=q)
  return(ff)

}


############################################################################
sieve.Pol <- function(X, J_n){
  return(unname(outer(X, 0:J_n, "^")));
}

############################################################################
#Trigonometric polynomials
#This function takes input vector x, and degree J_n and returns trigonometryic
#polynomials on [0,1] of degree J_n
sieve.TriPol <- function(X, J_n){
  if (J_n < 1){
    stop("J_n need to be at least 1");
  }
  basis <- rep(1,length(X));
  for(k in 1:J_n){
    basis <- cbind(basis, cos(2*k*pi*X), sin(2*k*pi*X));
  }
  return(unname(basis));
} #sieve.TriPol
#This function takes input vector x, degree J_n and order of derivative d and
#returns the derivative of the trigonometric polynomial basis
sieve.TriPol.D <- function(X, J_n, d){
  if (J_n < 1){
    stop("J_n need to be at least 1");
  }
  basis <- rep(0,length(X));
  if(d==1){
    for(k in 1:J_n){
      basis <- cbind(basis, -2*k*pi*sin(2*k*pi*X), 2*k*pi*cos(2*k*pi*X));
    } #end for loop
  } #end d==1
  else{
    for(k in 1:J_n){
      basis <- cbind(basis, -4*k^2*pi^2*cos(2*k*pi*X),
                     -4*k^2*pi^2*sin(2*k*pi*X));
    } #end for loop
  } #end d==1
  return(unname(basis));
} #end sieve.TriPol.D
############################################################################