構造方程式モデルと計量経済学

```{r}
library(dplyr)

n=100;G=8;K=7

Y=matrix(rpois(n*G,5),ncol=G)

Z=matrix(rpois(n*K,9),ncol=K)

P=solve(t(Z)%*%Z)%*%t(Z)%*%Y

H=t(Y-Z%*%P)%*%(Y-Z%*%P)



K1=2;G1=3;K2=K-K1

L=K2-G1

if(L>=0){

Z1=Z[,1:K1];Z2=Z[,-c(1:K1)]

Y1=Y[,1:G1];Y2=Y[,-c(1:G1)]

A11=t(Z1)%*%Z1

A12=t(Z1)%*%Z2

A21=t(Z2)%*%Z1

A22=t(Z2)%*%Z2

Z21=Z2-Z1%*%solve(t(Z1)%*%Z1)%*%(t(Z1)%*%Z2)

A=t(Z)%*%Z

A_star=rbind(t(Z1),t(Z21))%*%cbind(Z1,Z21)

P_star=solve(A_star)%*%rbind(t(Z1),t(Z21))%*%Y

H_star=t(Y)%*%Y-t(P_star)%*%A_star%*%P_star

g=t(Y)%*%(Z%*%solve(t(Z)%*%Z)%*%t(Z)-Z1%*%solve(t(Z1)%*%Z1)%*%t(Z1))%*%Y

ma1=t(cbind(Y,Z1))%*%Z%*%solve(t(Z)%*%Z)%*%t(Z)%*%cbind(Y,Z1)



f=function(param){
  
param=c(1,param)  
  
val1=as.numeric(t(param)%*%t(cbind(Y,Z1))%*%Z%*%solve(t(Z)%*%Z)%*%t(Z)%*%cbind(Y,Z1)%*%t(t(param)))  
  
val2=as.numeric(t(param)%*%t(cbind(Y,Z1))%*%cbind(Y,Z1)%*%t(t(param)))

return(val1/val2)

}


ite=20000

eta=0.01

h=0.01

params=rep(1,G+K1-1)

for(l in 1:ite){
  
vec=params

for(j in 1:length(vec)){
  
vec_sub=vec;vec_sub[j]=vec_sub[j]+h

params[j]=params[j]-eta*(f(vec_sub)-f(vec))/h
  
}
  
print(f(params))  

}

b2=-params[1:G1];r1=-params[-c(1:G1)]

}

```


```{r}

#p.45

library(dplyr)

n=100;G=8;K=7

Y=matrix(rpois(n*G,5),ncol=G)

Z=matrix(rpois(n*K,9),ncol=K)





K1=2;G1=3;K2=K-K1

L=K-G1-K1

Z1=t(Z[,1:K1]);

Y2=t(Y[,1:G1]);Y1=t(t(Y[,c(G1+1)]))

X=cbind(t(Y2),t(Z1))



#GMM推定量：P

func=function(s){

u=Y1-X%*%s

g=apply(Z*c(u),2,mean)

W_n=0

for(i in 1:n){
  
W_n=W_n+t(t(Z[i,]))%*%t(Z[i,])*u[i]^2  
  
}

W_n=W_n/n

return(as.numeric(t(g)%*%solve(W_n)%*%t(t(g))))

}


ite=10000

P=rep(1,ncol(X))

eta=0.01

h=0.01

for(k in 1:ite){
  
vec=P

for(j in 1:length(vec)){
  
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
  
P[j]=P[j]-eta*(func(vec_sub)-func(vec))/h

}

print(func(P))
  
}




mat1=array(0,dim=c(ncol(X),ncol(Z)))

mat2=array(0,dim=c(ncol(Z),1))

mat3=array(0,dim=c(ncol(Z),ncol(X)))

u=Y1-X%*%P;B1=array(0,dim=c(ncol(X),ncol(Z)));B2=array(0,dim=c(ncol(Z),ncol(X)))

W_n=array(0,dim=c(ncol(Z),ncol(Z)))

for(j in 1:n){

mat1=mat1+t(t(c(Z1[,j],Y2[,j])))%*%t(Z[j,])

mat2=mat2+t(t(Z[j,]))*Y1[j]

mat3=mat3+t(t(Z[j,]))%*%t(c(Z1[,j],Y2[,j]))

W_n=W_n+t(t(Z[i,]))%*%t(Z[i,])*u[i]^2

B1=B1+t(t(X[j,]))%*%t(Z[j,])

B2=B2+t(t(Z[j,]))%*%t(X[j,])

}

mat1=mat1/n;mat2=mat2/n;mat3=mat3/n;W_n=W_n/n;B1=B1/n;B2=B2/n

val1=mat1%*%(svd(W_n)$u%*%diag(1/svd(W_n)$d)%*%t(svd(W_n)$v))%*%mat2

val2=mat1%*%(svd(W_n)$u%*%diag(1/svd(W_n)$d)%*%t(svd(W_n)$v))%*%mat3

params=(svd(val2)$u%*%diag(1/svd(val2)$d)%*%t(svd(val2)$v))%*%val1

AGM=(B1%*%W_n%*%B2);AGM=svd(AGM)$u%*%diag(1/svd(AGM)$d)%*%t(svd(AGM)$v)

#H:delta=delta0

c=rep(1,length(params))

delta0=c*rep(0,length(params))

delta=c*params

sqrt(n)*(delta-delta0)/as.numeric(sqrt(t(c)%*%AGM%*%t(t(c))))

qnorm(1-0.05/2)

#H:R%*%P=c

R=diag(1,length(params))

c=rep(1,length(params))

svd(t(R)%*%AGM%*%R)$u%*%diag(1/svd(t(R)%*%AGM%*%R)$d)%*%t(svd(t(R)%*%AGM%*%R)$v)

f=n*t(R%*%params-c)%*%(svd(t(R)%*%AGM%*%R)$u%*%diag(1/svd(t(R)%*%AGM%*%R)$d)%*%t(svd(t(R)%*%AGM%*%R)$v))%*%(R%*%params-c)

qchisq(1-0.05,df=length(params))



```


```{r}

#p.51 推定方程式と経験尤度法

rm(list=ls())

library(dplyr)

n=30;G=8;K=7

Y=matrix(rpois(n*G,5),ncol=G)

Z=matrix(rpois(n*K,9),ncol=K)

K1=2;G1=3;K2=K-K1

L=K-G1-K1

Z1=t(Z[,1:K1]);

Y2=t(Y[,1:G1]);Y1=t(t(Y[,c(G1+1)]))

X=cbind(t(Y2),t(Z1))


log_lik=function(x){
  
p=x[1:n];x=x[-c(1:n)];

lam=x[1:ncol(Z)];x=x[-c(1:ncol(Z))]  

param=x;p=abs(p)/sum(abs(p));lam=abs(lam)

return(sum(lam*apply(Z*c(abs(Y1-X%*%param))*c(p),2,sum)))
  
}

params=rep(0.1,n+ncol(Z)+ncol(X))



ite=10000

h=0.01

eta=0.0001

log_lik_vec1=c()

log_lik_vec2=c()

for(l in 1:ite){
  
vec=params

for(j in 1:length(vec)){

vec_sub=vec;vec_sub[j]=vec_sub[j]+h  
    
params[j]=params[j]-eta*(log_lik(vec_sub)-log_lik(vec))/h  

}

params_sub=params

#print(log_lik(params))

p=params[1:n];params=params[-c(1:n)];

lam=params[1:ncol(Z)];params=params[-c(1:ncol(Z))]  

param=params;p=abs(p)/sum(abs(p));lam=abs(lam)

likli=-sum(log(1+lam*apply(Z*c((Y1-X%*%param))*c(p),2,sum)))
  
#print(paste0(likli,":",log_lik(params_sub)))

log_lik_vec1=c(log_lik_vec1,likli)

log_lik_vec2=c(log_lik_vec2,log_lik(params_sub))

params=params_sub

}



log_lik_vec1=log_lik_vec1[-c(1:(ite-100))]

log_lik_vec2=log_lik_vec2[-c(1:(ite-100))]

plot(log_lik_vec1)

```