LoginSignup
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@kozakai-ryouta(亮太 小酒井)

ベイジアン計量経済学入門　A.ゼルナー著　培風館

Last updated at Posted at 2022-01-09

y=c(0.699,0.32,-0.799,-0.927,0.373,-0.648,1.572,-0.319,2.049,-3.077)

mu_hat=mean(y)

n=length(y)

sig0=1;mu_a=6;sig_a=sqrt(200)

E_mu=(mu_hat*sig_a^2+mu_a*sig0/n)/(sig_a^2+sig0/n)

V_mu=(sig0*sig_a^2/n)/(sig_a^2+sig0/n)

hist(rnorm(10000,E_mu,V_mu))


# p.27

y=c(0.699,0.32,-0.799,-0.927,0.373,-0.648,1.572,-0.319,2.049,-3.077)

n=length(y)

mu_hat=mean(y)

vs=sum((y-mu_hat)^2)

p_sig_y=function(sigma){
  
  return(exp(-vs/(2*sigma^2))*sqrt(2*pi*sigma^2/n)*sigma^(-(n+1)))
  
}

plot(seq(0.01,5,0.01),p_sig_y(seq(0.01,5,0.01)))

sig_dat=data.frame(sig=seq(0.01,5,0.01),p=p_sig_y(seq(0.01,5,0.01)))

p_mu_y=function(mu){
  
  return((vs+n*(mu-mu_hat)^2)^(-(n+1)/2))
  
}

plot(seq(-5,5,0.01),p_mu_y(seq(-5,5,0.01)))

mu_dat=data.frame(mu=seq(-5,5,0.01),p=p_mu_y(seq(-5,5,0.01)))



# パレート分布の最尤推定 p.39

A=30

y=sample(A+1,200,50)

n=50

alpha=n/(sum(log(y))-n*log(A))

pareto=function(x){
  
return(n*log(x)+n*x*log(A)-(x+1)*sum(log(y)))
  
}

plot(seq(0.001,0.5,0.001),pareto(seq(0.001,0.5,0.001)))

G=prod(y)^(1/n)

a=log(G/A)

E_alpha=1/a

alpha_moment=function(r){
  
return(((a*n)^(-r))*gamma(n+r)/gamma(n))  
  
}

V_alpha=alpha_moment(2)-(alpha_moment(1))^2





# p.41

A=1

xt=c(10000,15000,25000,50000,100000,150000,500000)/10000

nt=c(171,222,160,220,48,138,45)

N=sum(nt)

log_G=sum(nt*log(xt))/N

as=log_G/A

func=function(x){
  
# return(exp(-as*N*x)*prod((1-(xt[1:(length(xt)-1)]/xt[2:length(xt)])^x)^nt[1:(length(xt)-1)])/x)
  
  
return((x*N*log(A)-x*sum(nt*log(xt))+sum(nt[1:(length(xt)-1)]*log((1-(xt[1:(length(xt)-1)]/xt[2:length(xt)])^x)))-log(x)))  
  
}

sequences=seq(0.001,2,0.001)

values=c()

for(j in 1:length(sequences)){
  
values=c(values,func(sequences[j]))  
  
}

values=values/L

plot(sequences,values)

# 勾配上昇法

alpha_hat=1

ite=10^4

eta=10^(-6)

h=0.001

for(l in 1:ite){
  
alpha_hat=alpha_hat+eta*(func(alpha_hat+h)-func(alpha_hat))/h

print(func(alpha_hat))
  
}




# p.66

data(anscombe)

X=anscombe$y4

Y=anscombe$y3

n=length(Y)

v=n-3

b2=sum((X-mean(X))*(Y-mean(Y)))/sum((X-mean(X))^2)

b1=mean(Y)-mean(X)*b2

s2=sum((Y-b1-b2*X)^2)/v

pb12=function(beta1,beta2){

value=(v*s2+n*(beta1-b1)^2+sum(X^2)*(beta2-b2)^2+2*(beta1-b1)*(beta2-b2)*sum(X))^(-n/2)

return(value)

}

pb1=function(beta1){

value=(v+sum((X-mean(X))^2)*(beta1-b1)^2)/(s2*sum(X^2)/n)

return(value^(-(v+1)/2))

}

pb2=function(beta2){

value=(v+sum((X-mean(X))^2)*(beta2-b2)^2)/s2

return(value^(-(v+1)/2))

}


b1+qt(1-0.05/2,df=v)*sqrt(sum((X-mean(X))^2)/(s2*sum(X^2)/n))^(-1)

b1-qt(1-0.05/2,df=v)*sqrt(sum((X-mean(X))^2)/(s2*sum(X^2)/n))^(-1)

b2+qt(1-0.05/2,df=v)*sqrt(sum((X-mean(X))^2))/s2

b2-qt(1-0.05/2,df=v)*sqrt(sum((X-mean(X))^2))/s2








# p.111 一部の係数が同一な場合の2本の回帰関係

n1=10
n2=20
ncol1=5
ncol2=6
ncol3=3

W1=matrix(sample(1:9,n1*ncol1,replace=T),nrow=n1)
Z1=matrix(sample(1:9,n2*ncol1,replace=T),nrow=n2)
W2=matrix(sample(1:9,n1*ncol2,replace=T),nrow=n1)
Z2=matrix(sample(20:30,n2*ncol3,replace=T),nrow=n2)

beta1=sample(1:9,ncol1,replace=T)
beta2=sample(1:9,ncol2,replace=T)
r2=sample(1:9,ncol3,replace=T)
Y1=W1%*%beta1+W2%*%beta2+rnorm(n1)
Y2=Z1%*%beta1+Z2%*%r2+rnorm(n2)
Y=c(Y1,Y2)

X1=cbind(W1,W2,array(0,dim=c(n1,ncol3)))
X2=cbind(Z1,array(0,dim=c(n2,ncol2)),Z2)
X=rbind(X1,X2)

sita=solve(t(X)%*%X)%*%t(X)%*%c(Y1,Y2)

v=n1+n2-ncol1-ncol2-ncol3
s2=sum((Y-X%*%sita)^2)/v

# betaの分布
pbeta=function(beta){
value=(v*s2+c(beta-sita)%*%(t(X)%*%X)%*%c(beta-sita))^(-(n1+n2)/2)
return(value)
}








# p.127 観測誤差モデル(EVM)
# 単回帰分析の従属変数・説明変数の双方の確率的な測定誤差を考慮に入れたモデル

n=100
b0=1
b=2
y1=5+rnorm(n,0,2)
y2=b0+b*5+rnorm(n,0,5)

ite=1000
sig1=4
sig2=2
guzai=rep(1,n)
eta=0.01

for(l in 1:ite){
  
guzai=(1-eta)*guzai+eta*((y1/sig1)+(b^2/sig2)*(y2-b0)/b)/(1/sig1+b^2/sig2)  
sita=sig1*b^2/sig2
w=(y2-b0)/b
sig1=(1-eta)*sig1+eta*sum((y1-w)^2)*sita^2/((1+sita^2)*n^2)  
sig2=(1-eta)*sig2+eta*sum((y1-w)^2)*b^2/(n*(1+sita^2))  

loglik=-n*log(sig1+1)-n*log(sig2+1)-sum((y1-guzai)^2)/(2*sig1)-sum((y2-b0-b*guzai)^2)/(2*sig2)

print(loglik)

}



# p.128

n=100
b0=1
b=2
y1=5+rnorm(n,0,2)
y2=b0+b*5+rnorm(n,0,5)
pthi=0.6

ite=10000
sig2=2
guzai=rep(1,n)
beta=1
beta0=2
eta=0.01

for(l in 1:ite){
beta0=beta0*(1-eta)+eta*sum((y2-beta0-beta*guzai))/(sig2^2)
beta=beta*(1-eta)+eta*sum((y2-beta0-beta*guzai)*guzai)/(sig2^2)  
guzai=(1-eta)*guzai+eta*(pthi*(y1-guzai)+beta*(y2-beta0-beta*guzai))/(sig2^2)
sig2=(1-eta)*sig2+eta*(-2*n/sig2+(pthi*sum((y1-guzai)^2)+sum((y2-beta0-beta*guzai)^2))/(sig2^3))

loglik=-2*n*log(sig2+1)-(pthi*sum((y1-guzai)^2)+sum((y2-beta0-beta*guzai)^2))/(2*sig2)

print(loglik)

}




# p.176

data(anscombe)

X=as.matrix(cbind(rep(1,nrow(anscombe)),anscombe[,7:8]))

V=anscombe[,5]

n=length(V)

beta=solve(t(X)%*%X)%*%t(X)%*%V

ite=10^3

for(l in 1:ite){

loglik=function(x){

z=log(V)+x*V

s=sum((z-X%*%beta)^2)/(n-3)

return(-n*log(s)/2+sum(log(1+x*V)))

}

sita=optim(1,loglik,control = list(fnscale = -1))$par

print(loglik(sita))

beta=solve(t(X)%*%X)%*%t(X)%*%(log(V)+sita*V)

}

r=exp(beta[1])

delta=1/(1+beta[2]/beta[3])

alpha=beta[3]/delta








# p.201 分布ラグモデル

# (7.48)

n=100
y0=1
x=rpois(n,40)
X=x[-c((length(x)-1):length(x))]
X1=x[-c(1,length(x))]
y=rpois(n,20)
Y1=y[-c(1,length(y))]
Y2=y[-c(1:2)]
Y=y[-c((length(y)-1):length(y))]

Z=cbind(Y1,X)
param=solve(t(Z)%*%Z)%*%t(Z)%*%Y
lam=param[1];a=param[2]
sig=sum((Y-Z%*%param)^2)/length(Y)


loglik=function(lambda,alpha,sigma){
G=diag(rep((1+lambda^2),length(Y)))
for(i in 1:(length(Y)-1)){
G[i,(i+1)]=-lambda
G[(i+1),i]=-lambda}
logliklihood=-log(det(G)+1)/2-length(Y)*log(sqrt(sig))-as.numeric(c(Y-Z%*%param)%*%solve(G)%*%c(Y-Z%*%param))/(2*sig)
return(logliklihood)}

ite=1000
eta=0.01
h=0.01

for(l in 1:ite){
G=diag(rep((1+lam^2),length(Y)))
for(i in 1:(length(Y)-1)){
G[i,(i+1)]=-lam
G[(i+1),i]=-lam}
lam=lam*(1-eta)+eta*(loglik(lam+h,a,sig)-loglik(lam,a,sig))/h
a=a*(1-eta)+eta*(loglik(lam,a+h,sig)-loglik(lam,a,sig))/h
sig=sig*(1-eta)+eta*as.numeric(c(Y-Z%*%param)%*%solve(G)%*%c(Y-Z%*%param))/length(Y)
print(loglik(lam,a,sig))
}






# p.201 分布ラグモデル

# (7.52)

n=100
y0=1
x=rpois(n,40)
X=x[-c((length(x)-1):length(x))]
X1=x[-c(1,length(x))]
y=rpois(n,20)
Y1=y[-c(1,length(y))]
Y2=y[-c(1:2)]
Y=y[-c((length(y)-1):length(y))]

Z=cbind(Y1,X)
sig=sum((Y-Z%*%param)^2)/length(Y)
param=solve(t(Z)%*%Z)%*%t(Z)%*%Y

ite=1000
eta=0.01

for(l in 1:ite){
param=param*(1-eta)+eta*solve(t(Z)%*%Z)%*%t(Z)%*%Y
loglik=-length(Y)*log(sqrt(sig))-sum((Y-Z%*%param)^2)/(2*sig)
sig=sig*(1-eta)+eta*sum((Y-Z%*%param)^2)/length(Y)
print(loglik)
}






# p.206

# 尤度p(lam,k,sig|C)の最大化

# Haavelmo

# 消費(投資)
C=c(39,60,42,52,47,51,45,60,39,41,22,17,27,33,48,51,33,46,54,100)

# 所得
Y=c(433,483,479,486,494,498,511,534,478,440,372,381,419,449,511,520,477,517,548,629)
# 年度
t=c(1922:1941)

ite=100

eta=0.001

lam=0.5
k=1
sig=1

y=Y[2:length(Y)]
c=C[2:length(C)]
cpast=C[-length(C)]

loglikvec=c()

for(l in 1:ite){
  
prek=k
prelam=lam
dlam=-sum((-cpast+k*y)*(c-lam*cpast-k*(1-lam)*y))/(sig^2)
dk=sum((y*(1-lam))*(c-lam*cpast-k*(1-lam)*y))/(sig^2)
sig=sqrt(sum((c-lam*cpast-k*(1-lam)*y)^2)/(length(t)+1))
lam=lam*(1-eta)+eta*dlam
k=k*(1-eta)+eta*dk
k=abs(k)/(abs(k)+abs(prek))
lam=abs(lam)/(abs(lam)+abs(prelam))
loglik=-(length(t)+1)*log(sig)-sum((c-lam*cpast-k*(1-lam)*y)^2)/(2*sig^2)
loglikvec=c(loglikvec,loglik)
print(loglik)
  
}





# p.206

# 尤度p(lam,k|C)の最大化(7.67)

# Haavelmo

# 消費(投資)
C=c(39,60,42,52,47,51,45,60,39,41,22,17,27,33,48,51,33,46,54,100)

# 所得
Y=c(433,483,479,486,494,498,511,534,478,440,372,381,419,449,511,520,477,517,548,629)
# 年度
t=c(1922:1941)

y=Y[2:length(Y)]
c=C[2:length(C)]
cpast=C[-length(C)]

# Gauss法
f=function(x){
return(x[1]*cpast+x[2]*(1-x[1])*y)
}
h=0.001
ite=10^2
sita=rep(0.1,2)
eta=10^(-2)
for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(length(c),length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub)-f(vec))/h
Z[,j]=Z_sub
}
beta=solve(t(Z)%*%Z)%*%t(Z)%*%(c-f(sita))
sita=c(sita+eta*beta)
print(sum((c-f(sita))^2))
}






# p.324 重回帰過程の単一期間制御

n1=10
n2=20
ncol1=5
ncol2=6
ncol3=3

W1=matrix(sample(1:9,n1*ncol1,replace=T),nrow=n1)
Z1=matrix(sample(1:9,n2*ncol1,replace=T),nrow=n2)
W2=matrix(sample(1:9,n1*ncol2,replace=T),nrow=n1)
Z2=matrix(sample(20:30,n2*ncol3,replace=T),nrow=n2)

beta1=sample(1:9,ncol1,replace=T)
beta2=sample(1:9,ncol2,replace=T)
r2=sample(1:9,ncol3,replace=T)
Y1=W1%*%beta1+W2%*%beta2+rnorm(n1)
Y2=Z1%*%beta1+Z2%*%r2+rnorm(n2)
Y=c(Y1,Y2)

X1=cbind(W1,W2,array(0,dim=c(n1,ncol3)))
X2=cbind(Z1,array(0,dim=c(n2,ncol2)),Z2)
X=rbind(X1,X2)

sita=solve(t(X)%*%X)%*%t(X)%*%c(Y1,Y2)

v=n1+n2-ncol1-ncol2-ncol3
s2=sum((Y-X%*%sita)^2)/v
sbar=v*s2/(v-2)

data(AirPassengers)
a=c(t(AirPassengers))

Wvec=c()
for(i in 1:length(a)){
w=t(X)%*%X%*%sita*a[i]/as.numeric(sbar+c(sita)%*%t(X)%*%X%*%c(sita))
Wvec=c(Wvec,w)
}

Wvec=matrix(Wvec,nrow=ncol1+ncol2+ncol3)
Zhat=t(Wvec)%*%sita


# 期待損失
Eza=c()
for(j in 1:length(a)){
Eza=c(Eza,sbar*(1+c(Wvec[,j])%*%solve(t(X)%*%X)%*%c(Wvec[,j]))+(a[j]-sum(Wvec[,j]*sita))^2)
}
0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?