---
title: "Factor IV Estimation in Conditional Moment Models with an application to Inflation Dynamics"
author: "Bertille Antoine and Xiaolin Sun"
date: "2024-03-08"
output:
  pdf_document: default
  html_document: default
  word_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


```{r}
library(MASS)
library(ivmodel)
library(sandwich)
```

 HAC Standard Errors
## define the functions 

```{r}
beta.smd.ts.het= function (nobs, X, Y, W, autocov){

W = data.frame(W)
  
b = dim(W)[2]
  
k = matrix(1,nrow = nobs,ncol = nobs)
  # Matrix of size n x length(x)
for (i in 1:b) {

  kx1 = sapply(W[,i], function(Xi) exp(-(W[,i] - Xi)^2 /2))
  k = k * kx1
}
 
 

ktil = k - k[1,1] * diag(nrow = nobs)

X = as.matrix(data.frame(X))

A = (1/nobs)*(1/(nobs))*t(X) %*% ktil %*% X
C = (1/nobs)*(1/(nobs))*t(X) %*% ktil %*% Y
beta =solve(A) %*% C

residual = Y -X %*% beta

ncolX = dim(X)[2]

mat1 = residual^2 %*% matrix(1,1,ncol = ncolX)
ktilX = ktil %*% X
B = t(( ktilX) * mat1) %*% ktilX

if (autocov==0) {
mat1 = mat1
B = B
  } else {
for (j in 1: autocov) {
 mat1.j = residual[(j+1):nobs] * residual[1:(nobs-j)] %*% matrix(1,1,ncol = ncolX)
 B.j = t(ktilX[(j+1):nobs,] * mat1.j) %*% ktilX[1:(nobs-j),]
 B =  B + (1 - j/(autocov+1)) * (B.j + t(B.j))
 }
}


varcov = solve(t(X)%*% ktil %*%  X)  %*% B %*% solve(t(X)%*% ktil %*%  X) 
var=diag(varcov)
se = sqrt(var)

seMat = matrix(se, dim(X)[2],1)
tstat = beta/seMat
  
output =data.frame(beta, seMat, tstat)

 return (output)
}
```


```{r}
beta.smd.band.ts.het= function (nobs, X, Y, W, c, autocov){

W = data.frame(W)
  
b = dim(W)[2]
  
k = matrix(1,nrow = nobs,ncol = nobs)
  # Matrix of size n x length(x)
for (i in 1:b) {

  kx1 = sapply(W[,i], function(Xi) exp(-(W[,i] - Xi)^2 /2))
  k = k * kx1
}
 
 
if (c == 0) {
ktil = k - k[1,1] * diag(nrow = nobs)
} else {

for (j in 1:c){
ktil = k - k[1,1] * diag(nrow = nobs)
  ktil[(1+j):nobs, 1:(nobs - j)] = ktil[(1+j):nobs, 1:(nobs - j)] - diag(diag(x = ktil[(1+j):nobs, 1:(nobs - j)]), nrow = (nobs - j))
}
}
ktil[upper.tri(ktil)] = t(ktil)[upper.tri(ktil)]

X = as.matrix(data.frame(X))

A = (1/nobs)*(1/(nobs))*t(X) %*% ktil %*% X
C = (1/nobs)*(1/(nobs))*t(X) %*% ktil %*% Y
beta =solve(A) %*% C

residual = Y -X %*% beta

ncolX = dim(X)[2]

mat1 = residual^2 %*% matrix(1,1,ncol = ncolX)
ktilX = ktil %*% X
B = t(( ktilX) * mat1) %*% ktilX

if (autocov==0) {
mat1 = mat1
B = B
  } else {
for (j in 1: autocov) {
 mat1.j = residual[(j+1):nobs] * residual[1:(nobs-j)] %*% matrix(1,1,ncol = ncolX)
 B.j = t(ktilX[(j+1):nobs,] * mat1.j) %*% ktilX[1:(nobs-j),]
 B =  B + (1 - j/(autocov+1)) * (B.j + t(B.j))
 }
}

varcov = solve(t(X)%*% ktil %*%  X)  %*% B %*% solve(t(X)%*% ktil %*%  X) 
var=diag(varcov)
se = sqrt(var)

seMat = matrix(se, dim(X)[2],1)
tstat = beta/seMat
  
output =data.frame(beta, seMat, tstat)

 return (output)
}
```


```{r}
beta.gmm.ts.eff.het= function (nobs, X, Z, Y, W, withZ, autocov){

 c= matrix(1,nobs,1)
  if (withZ==1) {
 Xmat = as.matrix(data.frame(c,X,Z))
 Zmat = as.matrix(data.frame(c,W,Z))
 } else {
 Xmat = as.matrix(data.frame(c,X))
 Zmat = as.matrix(data.frame(c,W))
 }
 
 
 proj = Zmat %*% solve(t(Zmat)%*% Zmat)  %*% t(Zmat)
 beta = solve(t(Xmat)%*% proj %*%  Xmat) %*% (t(Xmat)%*% proj %*% Y)
 residual = Y - Xmat %*% beta
 ncolX = dim(Zmat)[2] 
 mat1 = residual^2 %*% matrix(1,1,ncol = ncolX)
 A = t(Zmat * mat1) %*% Zmat
 
 beta.eff = solve(t(Xmat)%*% Zmat %*% solve(A)  %*% t(Zmat) %*%  Xmat) %*% (t(Xmat)%*% Zmat %*% solve(A) %*% t(Zmat) %*% Y)
 
residual = Y - Xmat %*% beta.eff
mat2 = (residual)^2 %*% matrix(1,1,ncol = ncolX)
A = t(Zmat * mat2) %*% Zmat
 
  if (autocov==0) {
mat2 = mat2
A = A
  } else {
 for (j in 1: autocov) {
 mat2.j = residual[(j+1):nobs] * residual[1:(nobs-j)] %*% matrix(1,1,ncol = ncolX)
 A.j = t(Zmat[(j+1):nobs,] * mat2.j) %*% Zmat[1:(nobs-j),]
 A =  A +  (1 - j/(autocov+1)) * (A.j+ t(A.j))
 }
}


 varcov = solve(t(Xmat)%*% Zmat %*% solve(A)  %*% t(Zmat) %*%  Xmat) 

 var=diag(varcov) #extract the diagonal
 se = sqrt(var)


 seMat = matrix(se, dim(Xmat)[2],1)
 tstat = beta/seMat
  
 output =data.frame(beta, seMat, tstat)
 
 return (output)
}
 
```


```{r}
gmm.jtest = function(beta, nobs, X, Z, Y, W, withZ){
  
  c= matrix(1,nobs,1)
  if (withZ==1) {
 Xmat = as.matrix(data.frame(c,X,Z))
 Zmat = as.matrix(data.frame(c,W,Z))
 } else {
 Xmat = as.matrix(data.frame(c,X))
 Zmat = as.matrix(data.frame(c,W))
 }

  beta.eff = beta
  ncolX = dim(Zmat)[2]
  residual = Y - Xmat %*% beta.eff
  mat1 = residual^2 %*% matrix(1,1,ncol = ncolX)
  A = t(Zmat * mat1) %*% Zmat
  
  Amat = t(Y - Xmat %*% beta.eff) %*% Zmat
  
  Jstat = Amat %*% solve(A) %*% t(Amat)
  return(Jstat)
}
```

```{r}
gmm.jtest.noconstant = function(beta, nobs, X, Z, Y, W, withZ){
  
  c= matrix(1,nobs,1)
  if (withZ==1) {
 Xmat = as.matrix(data.frame(X,Z))
 Zmat = as.matrix(data.frame(W,Z))
 } else {
 Xmat = as.matrix(data.frame(X))
 Zmat = as.matrix(data.frame(W))
 }

  beta.eff = beta
  ncolX = dim(Zmat)[2]
  residual = Y - Xmat %*% beta.eff
  mat1 = residual^2 %*% matrix(1,1,ncol = ncolX)
  A = t(Zmat * mat1) %*% Zmat
  
  Amat = t(Y - Xmat %*% beta.eff) %*% Zmat
  
  Jstat = Amat %*% solve(A) %*% t(Amat)
  return(Jstat)
}
```

```{r}
beta.Bgmm.eff.ts.het= function (nobs, X, Z, Y, W, withZ, Trf, autocov){
 Trf = round(Trf, digits = 0)  
b = dim(W)[2]

 
c_1 = Trf-1
c1 = matrix(1, nrow = c_1, ncol = 1)
 if (withZ==1) {
 Zmat1 = as.matrix(data.frame(c1,W[1:c_1,-b],Z[1:c_1,]))
} else {
 Zmat1 = as.matrix(data.frame(c1,W[1:c_1,-b]))
}

 c_2 = nobs-Trf+1
 c2 = matrix(1, nrow = c_2, ncol = 1)
 if (withZ==1) {
 Zmat2 = as.matrix(data.frame(c2,W[Trf:nobs, -b],Z[Trf:nobs, ]))
} else {
 Zmat2 = as.matrix(data.frame(c2,W[Trf:nobs, -b]))
}
 
 b = dim(Zmat1)[2]
 zeromat1 = matrix(0, nrow = c_1, ncol = b)

 zeromat2 = matrix(0, nrow = c_2, ncol = b)
 
 Zmat1 = cbind(Zmat1,zeromat1)
 Zmat2 = cbind(zeromat2,Zmat2)
 Zmat = rbind(Zmat1, Zmat2)
 
 
 c= matrix(1,nobs,1)
  if (withZ==1) {
 Xmat = as.matrix(data.frame(c,X,Z))
 } else {
 Xmat = as.matrix(data.frame(c,X))
 }
 
 proj = Zmat %*% solve(t(Zmat)%*% Zmat)  %*% t(Zmat) # it can be changed into proj = Zmat %*% t(Zmat)
 beta = solve(t(Xmat)%*% proj %*%  Xmat) %*% (t(Xmat)%*% proj %*% Y)
 residual = Y - Xmat %*% beta
 ncolX = dim(Zmat)[2] 
 mat1 = residual^2 %*% matrix(1,1,ncol = ncolX)
 A = t(Zmat * mat1) %*% Zmat
 
 beta.eff = solve(t(Xmat)%*% Zmat %*% solve(A)  %*% t(Zmat) %*%  Xmat) %*% (t(Xmat)%*% Zmat %*% solve(A) %*% t(Zmat) %*% Y)
 
residual = Y - Xmat %*% beta.eff
mat2 = (residual)^2 %*% matrix(1,1,ncol = ncolX)
A = t(Zmat * mat2) %*% Zmat
 
  if (autocov==0) {
mat2 = mat2
A = A
  } else {
 for (j in 1: autocov) {
 mat2.j = residual[(j+1):nobs] * residual[1:(nobs-j)] %*% matrix(1,1,ncol = ncolX)
 A.j = t(Zmat[(j+1):nobs,] * mat2.j) %*% Zmat[1:(nobs-j),]
 A =  A +  (1 - j/(autocov+1)) * (A.j+ t(A.j))
 }
}


 varcov = solve(t(Xmat)%*% Zmat %*% solve(A)  %*% t(Zmat) %*%  Xmat) 

 var=diag(varcov) #extract the diagonal
 se = sqrt(var)


 seMat = matrix(se, dim(Xmat)[2],1)
 tstat = beta/seMat
  
 output =data.frame(beta, seMat, tstat)
 
 return (output)
}
```
 

```{r}
beta.2sls.ts.het= function (nobs, X, Z, Y, W, withZ, autocov){
  
 c= matrix(1,nobs,1)
 if (withZ==1) {
 Xmat = as.matrix(data.frame(c,X,Z))
 Zmat = as.matrix(data.frame(c,W,Z))
 } else {
 Xmat = as.matrix(data.frame(c,X))
 Zmat = as.matrix(data.frame(c,W))
 }
 
 proj = Zmat %*% solve(t(Zmat)%*% Zmat)  %*% t(Zmat)
 beta = solve(t(Xmat)%*% proj %*%  Xmat) %*% (t(Xmat)%*% proj %*% Y)
 residual = Y - Xmat %*% beta
 
ncolX = dim(Zmat)[2]
sigma = mean(residual^2)

mat1 = residual^2 %*% matrix(1,1,ncol = ncolX)
A = t(Zmat * mat1) %*% Zmat

if (autocov==0) {
mat1 = mat1
A = A
  } else {
 for (j in 1: autocov) {
 mat1.j = residual[(j+1):nobs] * residual[1:(nobs-j)] %*% matrix(1,1,ncol = ncolX)
 A.j = t(Zmat[(j+1):nobs,] * mat1.j) %*% Zmat[1:(nobs-j),]
 A =  A +  (1 - j/(autocov+1)) * (A.j+ t(A.j))
 }
}

varcov = solve(t(Xmat) %*% proj %*%  Xmat) %*% t(Xmat) %*% Zmat %*% solve(t(Zmat) %*% Zmat)  %*% A %*% solve(t(Zmat) %*% Zmat)  %*% t(Zmat)  %*%  Xmat %*% solve(t(Xmat) %*% proj %*%  Xmat) 
 


var=diag(varcov) #extract the diagonal
se = sqrt(var)

seMat = matrix(se, dim(Xmat)[2],1)
tstat = beta/seMat
  
output =data.frame(beta, seMat, tstat)
 
 return (output)
}
```



## example



```{r}

mu_2 = c(0,0) 



beta = 0


cov_uv = 0.6
sigma2 = matrix(c(1,cov_uv,cov_uv, 1),2,2)
sigma2

n_ins = 50 
n_loadings = n_ins*2
set.seed(1)
a = rnorm(n_loadings, 1,1)
set.seed(NULL)
mu = matrix(a, nrow = n_ins, ncol = 2)
mu_3 = rep(0, n_ins)
sigma3 = diag(n_ins)

nobsc = c(200, 2000, 200, 2000)
repe = 5000
probz2 = c(0.2, 0.2, 0.05, 0.05)



m=3


nobs = nobsc[m]
lambda = probz2[m]
T_break = (1- lambda)*nobs+1
set.seed(2024)




vandu = mvrnorm(nobs,mu_2,sigma2)

F1 =  runif(nobs,min = -2, max = 2)


F2 = matrix(data=0,nrow = nobs,ncol = 1)
F2[T_break: nobs,1] = 1

Factor = as.matrix(data.frame(F1,F2))
error = mvrnorm(nobs,mu_3,sigma3)

Instruments = Factor %*% t(mu) + error

Y = 10 *(2*F2-1)* (F1 - 2*(F1^3)/5) + vandu[,1]

sig1 = matrix(NA, nobs, 1)
sig1[1] = abs(rnorm(1, mean = 0, sd = 1))
 for (l in 2:nobs) {
  sig1[l] = sqrt(0.1+ 0.6* (sig1[l-1]^2) * (vandu[l-1,2]^2) + 0.3 * (sig1[l-1]^2))
 }


y = Y * beta + vandu[,2]* sig1 

Xvar = Y

pr.ins= prcomp(Instruments,scale = T)

pc1 = pr.ins$x[,1]
pc2 = pr.ins$x[,2]
pc3 = pr.ins$x[,3]

W = as.matrix(data.frame(pc1,pc2,pc3))



 #nobs: number of observations
  #X: the endogenous variable. Xvar is a variable or a vector of variables
  #Y: outcome. y
  #W: the estimated factors.
  #c: the width of the band that eliminates pairs of observations (say s and t) that are not only equal to each other, but also too close to each other. Defined in Equation (30) in the Supplementary Appendix
  #autocov: the truncation parameter in HAC Standard Errors to be chosen. A rule of thumb for choosing m is m = ceiling(0.75*nobs^(1/3))-1. 

m = ceiling(0.75*nobs^(1/3))-1
c = matrix(1,nrow = nobs,ncol = 1)
#### HAC Standard Errors for all estimators
beta.smd.ts.het(nobs,X=Xvar,Y=y,W=W[,1],autocov = m) #FSMD with HAC Standard Errors
beta.smd.ts.het(nobs,X=as.matrix(data.frame(c,Xvar)),Y=y,W=W[,1],autocov = m)
beta.smd.band.ts.het(nobs, X=Xvar,Y=y,W=W[,1],c=1, autocov = m)
beta.smd.band.ts.het(nobs, X = as.matrix(data.frame(c,Xvar)), Y = y, W=W[,1], c=1, autocov = m) 

beta.gmm.ts.eff.het(nobs, X = Xvar, Z=Xvar, Y = y, W=W[,1], withZ=0,autocov = m)#GMM with HAC Standard Errors
beta.2sls.ts.het(nobs, X = Xvar, Z=Xvar, Y = y, W=W[,1], withZ=0,autocov = m)

beta.smd.ts.het(nobs,X=Xvar,Y=y,W=W[,1:2],autocov = 1) #FSMD with HAC Standard Errors
beta.smd.ts.het(nobs,X=as.matrix(data.frame(c,Xvar)),Y=y,W=W[,1:2],autocov = m)
beta.smd.band.ts.het(nobs, X=Xvar,Y=y,W=W[,1:2],c=1, autocov = m)
beta.smd.band.ts.het(nobs, X = as.matrix(data.frame(c,Xvar)), Y = y, W=W[,1:2], c=1, autocov = m) 

beta.gmm.ts.eff.het(nobs, X = Xvar, Z=Xvar, Y = y, W=W[,1:2], withZ=0,autocov = m)#GMM with HAC Standard Errors
beta.2sls.ts.het(nobs, X = Xvar, Z=Xvar, Y = y, W=W[,1:2], withZ=0,autocov = m)

beta.smd.ts.het(nobs,X=Xvar,Y=y,W=W,autocov = 1) #FSMD with HAC Standard Errors
beta.smd.ts.het(nobs,X=as.matrix(data.frame(c,Xvar)),Y=y,W=W,autocov = m)
beta.smd.band.ts.het(nobs, X=Xvar,Y=y,W=W,c=1, autocov = m)
beta.smd.band.ts.het(nobs, X = as.matrix(data.frame(c,Xvar)), Y = y, W=W, c=1, autocov = m) 


beta.gmm.ts.eff.het(nobs, X = Xvar, Z=Xvar, Y = y, W=W, withZ=0,autocov = m)#GMM with HAC Standard Errors
beta.2sls.ts.het(nobs, X = Xvar, Z=Xvar, Y = y, W=W, withZ=0,autocov = m)

```