# 動学的パネルデータ分析
# p.46
library(dplyr)
N=50;t=12
data=data.frame(no=1:144,data.frame(AirPassengers))
D_t=array(0,dim=c(t-1,t))
for(j in 1:nrow(D_t)){
D_t[j,(j:(j+1))]=c(-1,1)
}
Q_t=t(D_t)%*%solve(D_t%*%t(D_t))%*%D_t
# Q_t=diag(1,t)-t(t(rep(1,t)))%*%t(rep(1,t))/t
up=0;under=0
upper=0;und=0
for(j in 1:(N)){
yi=c(data$AirPassengers[(j+1):(j+t)])
yi_past=c(data$AirPassengers[j:(j+t-1)])
up=up+as.numeric(t(yi_past)%*%Q_t%*%t(t(yi)))
under=under+as.numeric(t(yi_past)%*%Q_t%*%t(t(yi_past)))
upper=upper+sum(yi*yi_past)
und=und+sum(yi_past*yi_past)
}
alpha_LSDV=up/under
alpha_POLS=upper/und
alpha_HK=(t+1)*alpha_LSDV/t+1/t
library(dplyr)
N1=90;N2=50;t=50
data=data.frame(no=1:144,data.frame(AirPassengers))
alpha_HK_vec=c();alpha_IV_D_vec=c();alpha_POLS_vec=c()
for(l in N2:N1){
D_t=array(0,dim=c(t-1,t))
for(j in 1:nrow(D_t)){
D_t[j,(j:(j+1))]=c(-1,1)
}
Q_t=t(D_t)%*%solve(D_t%*%t(D_t))%*%D_t
# Q_t=diag(1,t)-t(t(rep(1,t)))%*%t(rep(1,t))/t
up=0;under=0
upper=0;und=0
up_IV=0;under_IV=0
up_IV_BO=0;under_IV_BO=0
for(j in 1:(l)){
yi=c(data$AirPassengers[(j+1):(j+t)])
yi_past=c(data$AirPassengers[j:(j+t-1)])
up=up+as.numeric(t(yi_past)%*%Q_t%*%t(t(yi)))
under=under+as.numeric(t(yi_past)%*%Q_t%*%t(t(yi_past)))
upper=upper+sum(yi*yi_past)
und=und+sum(yi_past*yi_past)
if(j>2){
yi_pp=c(data$AirPassengers[(j-1):(j+t-2)])
dyi=yi-yi_past
yi2=c(data$AirPassengers[(j):(j+t-1)])
yi_past2=c(data$AirPassengers[(j-1):(j+t-2)])
dyi2=yi2-yi_past2
up_IV=up_IV+sum(yi_pp*dyi)
under_IV=under_IV+sum(dyi2*yi_pp)
}
}
alpha_LSDV=up/under
alpha_POLS=upper/und
alpha_HK=(t+1)*alpha_LSDV/t+1/t
alpha_IV_BO=
alpha_IV_D=up_IV/under_IV
alpha_POLS_vec=c(alpha_POLS_vec,alpha_POLS)
alpha_IV_D_vec=c(alpha_IV_D_vec,up_IV/under_IV)
alpha_IV_BO_vec=c(alpha_IV_BO_vec,)
alpha_HK_vec=c(alpha_HK_vec,alpha_HK)
}
plot(alpha_HK_vec)
alpha=mean(alpha_HK_vec)
cost=function(a){
return((var(sqrt(N*t)*(alpha_HK_vec-a))-(1-a^2))^2+(a-mean(alpha_HK_vec))^2 )
}
ite=100000
eta=0.0001;h=0.01
for(l in 1:ite){
alpha=alpha-eta*(cost(alpha+h)-cost(alpha))/h
print(cost(alpha))
}
# p.60
# GMM推定量
library(dplyr)
t=12;N=50
data=data.frame(no=1:144,data.frame(AirPassengers))
D_t=array(0,dim=c(t-2,t-1))
for(j in 1:(nrow(D_t))){
D_t[j,(j:(j+1))]=c(-1,1)
}
V_D_L2=array(0,dim=c((t-1)*(t-2)/2,(t-1)*(t-2)/2))
past_left=array(0,dim=c(1,(t-1)*(t-2)/2))
past_right=array(0,dim=c((t-1)*(t-2)/2,1))
right=array(0,dim=c((t-1)*(t-2)/2,1))
Z2=array(0,dim=c((t-1)*(t-2)/2,(t-1)*(t-2)/2))
for(j in 1:N){
y_vec=c(data$AirPassengers[j:(j+t-1)])
mat=array(0,dim=c((t-2),1));mat[1,1]=y_vec[1]
for(i in 2:(t-2)){
Z_i_sub=array(0,dim=c(t-2,i))
Z_i_sub[i,]=c(y_vec[1:i])
mat=cbind(mat,Z_i_sub)
}
Z_i_L2=t(mat)
Z2=Z2+Z_i_L2%*%t(Z_i_L2)
V_D_L2=V_D_L2+Z_i_L2%*%D_t%*%t(D_t)%*%t(Z_i_L2)
if(j>2){
y_past_vec=c(data$AirPassengers[(j-2):(j-2+t-3)])
y_vec=c(data$AirPassengers[(j-1):(j+t-4)])
y_past_vec2=c(data$AirPassengers[(j-1):(j-2+t-2)])
y_vec2=c(data$AirPassengers[j:(j+t-3)])
dyi=y_vec-y_past_vec;dyi_present=y_vec2-y_past_vec2
past_left=past_left+t(dyi)%*%t(Z_i_L2)
past_right=past_right+(Z_i_L2)%*%t(t(dyi))
right=right+(Z_i_L2)%*%t(t(dyi_present))
}
}
alpha_GMM=as.numeric(past_left%*%V_D_L2%*%right)/as.numeric(past_left%*%V_D_L2%*%past_right)
alpha_GMM_F=as.numeric(past_left%*%solve(Z2)%*%right)/as.numeric(past_left%*%solve(Z2)%*%past_right)