-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathGMM_Est.R
97 lines (80 loc) · 3.1 KB
/
GMM_Est.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
GMM_Est = function(Y,X,Z,k,l,W,r=0){
# require tidyverse
#X : instrument, Z: regressors, W: weight, r : restriction on coeff(1, same coeff)
#k: seq of number of instruments, l:seq of number of regressors
if(r==1){ # same coeff
nrowY = nrow(Y) #obs
ncolY = ncol(Y) #dim
nrowX = nrow(X)
ncolX = ncol(X) # sum(K)
nrowZ = nrowX
ncolZ = ncol(Z) # sum(L)
for(i in 1:length(k)){
if(i==1){
block = X[,1:k[i]]
block01 = Z[,1:l[i]]
EsigmaXY = 1/nrowX*t(block)%*%Y[1:nrowX,]
ESSigmaXZ = 1/nrowX*t(block)%*%block01
}else{
block = X[,(sum(k[1:(i-1)])+1):sum(k[1:i])]
block01 = Z[,(sum(l[1:(i-1)])+1):sum(l[1:i])]
EsigmaXY=rbind(EsigmaXY,1/nrowX*t(block)%*%Y[(nrowX*(i-1)+1):(nrowX*i),])
ESSigmaXZ = rbind(ESSigmaXZ,1/nrowX*t(block)%*%block01)
}
}
GMM.coeff = solve(t(ESSigmaXZ)%*%W%*%ESSigmaXZ)%*%t(ESSigmaXZ)%*%W%*%EsigmaXY
for(i in 1:length(k)){ # i번째 식, # k[i]: i번째 식의 차원
if(i==1){
Error = Y[1:nrowZ,] - Z[,1:l[i]]%*%GMM.coeff #첫번째 식의 에러
g = 1/nrowZ*diag(as.numeric(Error))%*%(X[,1:k[i]])
}else{
ErrorNext = Y[(nrowZ*(i-1)+1):(nrowZ*i),] -
Z[,(sum(l[1:(i-1)])+1):sum(l[1:i])]%*%GMM.coeff
gNext = 1/nrowZ*diag(as.numeric(ErrorNext))%*%(X[,(sum(k[1:(i-1)])+1):sum(k[1:i])])
Error = cbind(Error,ErrorNext)
g = cbind(g,gNext)
}
}
W01 = solve(1/nrowZ*t(g)%*%g)
return(list(GMM.coeff,W01))
#return(GMM.coeff)
}
else{
nrowY = nrow(Y) #obs
ncolY = ncol(Y) #dim
nrowX = nrow(X)
ncolX = ncol(X) # K
#nrowZ = nrow(Z)
ncolZ = ncol(Z) # L
for(i in 1:length(k)){
if(i==1){
block = X[,1:k[i]]
block01 = Z[,1:l[i]]
EsigmaXY = 1/nrowX*t(block)%*%Y[1:nrowX,]
#ESSigmaXZ = 1/nrowX*t(block)%*%block01
ESSigmaXZ = diag(ESSigmaXZ,1/nrowX*t(block)%*%block01)
}else{
block = X[,(sum(k[1:(i-1)])+1):sum(k[1:i])]
block01 = Z[,(sum(l[1:(i-1)])+1):sum(l[1:i])]
EsigmaXY=rbind(EsigmaXY,1/nrowX*t(block)%*%Y[(nrowX*(i-1)+1):(nrowX*i),])
#ESSigmaXZ = rbind(ESSigmaXY,1/nrowX*t(block)%*%block01)
ESSigmaXZ = bdiag(ESSigmaXZ,1/nrowX*t(block)%*%block01)
}
}
GMM.coeff = as.matrix(solve(t(ESSigmaXZ)%*%W%*%ESSigmaXZ )%*%t(ESSigmaXZ)%*%W %*%EsigmaXY)
for(i in 1:length(k)){ # i번째 식, # k[i]: i번째 식의 차원
if(i==1){
Error = Y[1:nrowZ,] - Z[,1:l[i]]%*%GMM.coeff[1:l[i],1] #첫번째 식의 에러
g = 1/nrowZ*diag(as.numeric(Error))%*%(X[,1:k[i]])
}else{
ErrorNext = Y[(nrowZ*(i-1)+1):(nrowZ*i),] -
Z[,(sum(l[1:(i-1)])+1):sum(l[1:i])]%*%GMM.coeff[(sum(l[1:(i-1)])+1):sum(l[1:i]),1]
gNext = 1/nrowZ*diag(as.numeric(ErrorNext))%*%(X[,(sum(k[1:(i-1)])+1):sum(k[1:i])])
Error = cbind(Error,ErrorNext)
g = cbind(g,gNext)
}
}
W01 = solve(1/nrowZ*t(g)%*%g)
return(GMM.coeff,W01)
}
}