固定効果モデルサンプル（動学的パネルデータ分析 千木良 弘朗 , 早川和彦, 山本拓先生 知泉書館）

#動学的パネルデータ分析　p.9~

 library(plm);library(dplyr)

data("EmplUK",package="plm")

N=max(sort(unique(EmplUK$firm)))

EmplUK$year=EmplUK$year-1976

x_ave=EmplUK %>% group_by(firm) %>% summarise(sector=mean(sector),emp=mean(emp),wage=mean(wage),capital=mean(capital),output=mean(output))

y_ave=x_ave %>% select(firm,output)

x_ave=x_ave %>% select(-output)

b_vec=c();etas=c();sigma_vec=c()

result_data=data.frame(firm=0,Y_Q=0,predict_Q=0,Y=0,predict=0)

for(j in 1:N){

X_sub=EmplUK %>% filter(firm==j)  

X_mat=as.matrix(X_sub[,colnames(X_sub) %in% c("emp","wage","capital")])

Y_sub=X_sub$output

t=nrow(X_mat)

D_t=array(0,dim=c(max(t)-1,max(t)))

for(k in 1:nrow(D_t)){

D_t[k,k:(k+1)]=c(1,-1)  

}

Q_t=t(D_t)%*%solve(D_t%*%t(D_t))%*%D_t

diff_x=t(X_mat)%*%Q_t%*%X_mat

diff_xy=t(X_mat)%*%Q_t%*%Y_sub

b=solve(diff_x)%*%diff_xy

eta=sum(Y_sub-X_mat%*%b)/t

b_vec=c(b_vec,b);etas=c(etas,eta)

res=data.frame(firm=j,Y_Q=Q_t%*%Y_sub,predict_Q=Q_t%*%X_mat%*%b,Y=Y_sub,predict=X_mat%*%b+eta)

result_data=rbind(result_data,res)

if(j==1){

diff_x_vec=diff_x  

diff_xy_vec=diff_xy  

}else{


diff_x_vec=diff_x_vec+diff_x

diff_xy_vec=diff_xy_vec+diff_xy

}

}


beta=solve(diff_x_vec)%*%diff_xy_vec

eta_hat=c()

result=data.frame(firm=0,Y_Q=0,predict_Q=0,Y=0,predict=0)

for(j in 1:N){

X_sub=EmplUK %>% filter(firm==j)  

X_mat=as.matrix(X_sub[,colnames(X_sub) %in% c("emp","wage","capital")])

Y_sub=X_sub$output

t=nrow(X_mat)

D_t=array(0,dim=c(max(t)-1,max(t)))

for(k in 1:nrow(D_t)){

D_t[k,k:(k+1)]=c(1,-1)  

}

eta=sum((Y_sub)-X_mat%*%beta)/t

eta_hat=c(eta_hat,eta)

Q_t=t(D_t)%*%solve(D_t%*%t(D_t))%*%D_t

res=data.frame(firm=j,Y_Q=Q_t%*%Y_sub,predict_Q=Q_t%*%X_mat%*%beta,Y=Y_sub,predict=X_mat%*%beta+eta)


result=rbind(result,res)

if(j==1){

sigma=t(Q_t%*%Y_sub-Q_t%*%X_mat%*%beta)%*%(Q_t%*%Y_sub-Q_t%*%X_mat%*%beta)  

}else{


sigma=sigma+t(Q_t%*%Y_sub-Q_t%*%X_mat%*%beta)%*%(Q_t%*%Y_sub-Q_t%*%X_mat%*%beta)  


}

}

sigma=sigma/(N*(t-1)-3)

result=tail(result,nrow(result)-1)

result_data=tail(result_data,nrow(result_data)-1)

plot(c(1:nrow(result)),result$Y,type="o",col=2,xlim=c(1,nrow(result)),ylim=c(min(c(result$Y,result$predict)),max(c(result$Y,result$predict))),xlab="index",ylab="values")


par(new=T)

plot(c(1:nrow(result)),result$predict,type="o",col=3,xlim=c(1,nrow(result)),ylim=c(min(c(result$Y,result$predict)),max(c(result$Y,result$predict))),xlab="index",ylab="values")

plot(c(1:nrow(result_data)),result_data$Y,type="o",col=2,xlim=c(1,nrow(result_data)),ylim=c(min(c(result_data$Y,result_data$predict)),max(c(result_data$Y,result_data$predict))),xlab="index",ylab="values")


par(new=T)

plot(c(1:nrow(result_data)),result_data$predict,type="o",col=3,xlim=c(1,nrow(result_data)),ylim=c(min(c(result_data$Y,result_data$predict)),max(c(result_data$Y,result_data$predict))),xlab="index",ylab="values")

 library(plm);library(dplyr)

data("EmplUK",package="plm")

N=max(sort(unique(EmplUK$firm)))

EmplUK$year=EmplUK$year-1976

x_ave=EmplUK %>% group_by(firm) %>% summarise(sector=mean(sector),emp=mean(emp),wage=mean(wage),capital=mean(capital),output=mean(output))

y_ave=x_ave %>% select(firm,output)

x_ave=x_ave %>% select(-output)

b_vec=c();etas=c();sigma_vec=c()

result_data=data.frame(firm=0,Y_Q=0,predict_Q=0,Y=0,predict=0)

for(j in 1:N){

X_sub=EmplUK %>% filter(firm==j)  

X_mat=as.matrix(X_sub[,colnames(X_sub) %in% c("emp","wage","capital")])

Y_sub=X_sub$output

t=nrow(X_mat)

D_t=array(0,dim=c(max(t)-1,max(t)))

for(k in 1:nrow(D_t)){

D_t[k,k:(k+1)]=c(1,-1)  

}

Q_t=t(D_t)%*%solve(D_t%*%t(D_t))%*%D_t

diff_x=t(X_mat)%*%Q_t%*%X_mat

diff_xy=t(X_mat)%*%Q_t%*%Y_sub

b=solve(diff_x)%*%diff_xy

eta=sum(Y_sub-X_mat%*%b)/t

b_vec=c(b_vec,b);etas=c(etas,eta)

res=data.frame(firm=j,Y_Q=Q_t%*%Y_sub,predict_Q=Q_t%*%X_mat%*%b,Y=Y_sub,predict=X_mat%*%b+eta)

result_data=rbind(result_data,res)

if(j==1){

diff_x_vec=diff_x  

diff_xy_vec=diff_xy  

}else{


diff_x_vec=diff_x_vec+diff_x

diff_xy_vec=diff_xy_vec+diff_xy

}

}


b_vec=matrix(b_vec,nrow=ncol(X_mat))

beta=solve(diff_x_vec)%*%diff_xy_vec


eta_hat=c()

result=data.frame(firm=0,Y_Q=0,predict_Q=0,Y=0,predict=0)

sigma_sum=c()

for(j in 1:N){

X_sub=EmplUK %>% filter(firm==j)  

X_mat=as.matrix(X_sub[,colnames(X_sub) %in% c("emp","wage","capital")])

Y_sub=X_sub$output

t=nrow(X_mat)

D_t=array(0,dim=c(max(t)-1,max(t)))

for(k in 1:nrow(D_t)){

D_t[k,k:(k+1)]=c(1,-1)  

}

sigma_sum=c(sigma_sum, t(Y_sub-X_mat%*%beta)%*%(Y_sub-X_mat%*%beta)/(N*t-ncol(X_mat)))

eta=sum((Y_sub)-X_mat%*%beta)/t

eta_hat=c(eta_hat,eta)

Q_t=t(D_t)%*%solve(D_t%*%t(D_t))%*%D_t

res=data.frame(firm=j,Y_Q=Q_t%*%Y_sub,predict_Q=Q_t%*%X_mat%*%beta,Y=Y_sub,predict=X_mat%*%beta+eta)


result=rbind(result,res)

if(j==1){

sigma=t(Q_t%*%Y_sub-Q_t%*%X_mat%*%beta)%*%(Q_t%*%Y_sub-Q_t%*%X_mat%*%beta)  

}else{


sigma=sigma+t(Q_t%*%Y_sub-Q_t%*%X_mat%*%beta)%*%(Q_t%*%Y_sub-Q_t%*%X_mat%*%beta)  


}

}

sigma=sigma/(N*(t-1)-3)

#固定効果での回帰係数の推定結果

result=tail(result,nrow(result)-1)

#固定効果での各firmごとの回帰係数の推定結果

result_data=tail(result_data,nrow(result_data)-1)



sigma_eta=sum(sigma_sum)-sigma

sigma_v=sigma

#変量効果の回帰係数の推定

for(j in 1:N){

X_sub=EmplUK %>% filter(firm==j)  

X_mat=as.matrix(X_sub[,colnames(X_sub) %in% c("emp","wage","capital")])

Y_sub=X_sub$output

t=nrow(X_mat)

omega_t=array(sigma_eta,dim=c(t,t))+diag(as.numeric(sigma_v),t)

if(j==1){

diff_xx=t(X_mat)%*%solve(omega_t)%*%X_mat

diff_xy=t(X_mat)%*%solve(omega_t)%*%Y_sub

}else{

diff_xx=diff_xx+t(X_mat)%*%solve(omega_t)%*%X_mat

diff_xy=diff_xy+t(X_mat)%*%solve(omega_t)%*%Y_sub


}

}

beta_re=solve(diff_xx)%*%diff_xy

V_beta_re=solve(diff_xx)



result_random=data.frame(firm=0,Y=0,predict=0)

r_vec=c();t_vec=c()

for(j in 1:N){

X_sub=EmplUK %>% filter(firm==j)  

X_mat=as.matrix(X_sub[,colnames(X_sub) %in% c("emp","wage","capital")])

Z=array(0,dim=c(nrow(X_mat),ncol(X_mat)))

for(i in 1:nrow(X_mat)){

Z[i,sample(c(1:ncol(X_mat)),1)]=1  

}

Y_sub=X_sub$output

diff=t(t(Y_sub))-as.matrix(X_mat)%*%beta_re

t=nrow(X_mat)

r=c(t(Z)%*%diff)/apply(Z,2,sum)

r[is.nan(r)>0]=0

r_vec=c(r_vec,r)

res=data.frame(firm=j,Y=Y_sub,predict=X_mat%*%beta_re+Z%*%r)

result_random=rbind(result_random,res)

t_vec=c(t_vec,t)

if(j==1){

Z_mat=Z  



}else{

Z_mat=rbind(Z_mat,Z)  

}

}

result_random=tail(result_random,nrow(result_random)-1)

cor(result_random$Y,result_random$predict)

plot(c(1:nrow(result_random)),result_random$Y,type="o",col=2,xlim=c(1,nrow(result_random)),ylim=c(min(c(result_random$Y,result_random$predict)),max(c(result_random$Y,result_random$predict))),xlab="index",ylab="values")


par(new=T)

plot(c(1:nrow(result_random)),result_random$predict,type="o",col=3,xlim=c(1,nrow(result_random)),ylim=c(min(c(result_random$Y,result_random$predict)),max(c(result_random$Y,result_random$predict))),xlab="index",ylab="values")