Go to Qiita Advent Calendar 2024 Top
LoginSignup
2
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(亮太 小酒井)

計量経済学の理論と応用　蓑谷千凰彦先生　日本評論社

Last updated at Posted at 2021-05-21


#p.103 不均一分散の場合

#勾配上昇法

data("AirPassengers")

p=3

data=data.frame(no=1:length(AirPassengers),values=c(AirPassengers))

r=4;dim=1

values=data$values

X=array(0,dim=c((nrow(data)-r),dim*r))

Y=c()

for(j in 1:(nrow(data)-r)){

vec=c()

for(i in 1:dim){

vec=c(vec,values[j:(j+r-1)]^i)  

}

X[j,]=vec

Y=c(Y,(values[r+j]))

}

beta=rep(0.001,ncol(X))

delta=rep(1,2)*10^(-5)

w=cbind(rep(1,nrow(X)),X%*%beta)

ite=10^6

eta=10^(-8)

for(l in 1:ite){

dbeta=t(X)%*%c(exp(-w%*%delta)*(Y-X%*%beta))

ddelta=-apply(w,2,sum)/2+t(w)%*%c(exp(-w%*%delta)*(Y-X%*%beta)^2)/2

beta=beta+eta*dbeta

delta=delta+eta*ddelta

loglik=-sum(w%*%delta)/2-sum(exp(-w%*%delta)*(Y-X%*%beta)^2)/2

print(loglik)

w=cbind(rep(1,nrow(X)),X%*%beta)

}

vt=(Y-X%*%beta)^2/exp(delta[1])-1

w_mat=array(0,dim=c(length(delta),length(delta)))

for(j in 1:nrow(w)){

w_mat=w_mat+t(t(c(w[j,])))%*%t(c(w[j,]))

}

nrow(X)*t(t(w)%*%vt)%*%solve(w_mat)%*%(t(w)%*%vt)/sum(vt^2)

qchisq(1-0.05,df=1)





#p.103 不均一分散の場合

#勾配上昇法 （予測子がexp(X%*%beta)の場合）

data("AirPassengers")

p=3

data=data.frame(no=1:length(AirPassengers),values=c(AirPassengers))

r=4;dim=1

values=data$values/100

X=array(0,dim=c((nrow(data)-r),dim*r))

Y=c()

for(j in 1:(nrow(data)-r)){

vec=c()

for(i in 1:dim){

vec=c(vec,values[j:(j+r-1)]^i)  

}

X[j,]=vec

Y=c(Y,(values[r+j]))

}



beta=rep(0.001,ncol(X))

delta=rep(1,2)*10^(-2)

w=cbind(rep(1,nrow(X)),X%*%beta)

loglik=-sum(w%*%delta)/2-sum(exp(-w%*%delta)*(Y-exp(X%*%beta))^2)/2

ite=10^7

eta=10^(-7)

for(l in 1:ite){

dbeta=t(X)%*%c(exp(-w%*%delta)*(Y-exp(X%*%beta))*exp(X%*%beta))

ddelta=-apply(w,2,sum)/2+t(w)%*%c(exp(-w%*%delta)*(Y-exp(X%*%beta))^2)/2

beta=beta+eta*dbeta

delta=delta+eta*ddelta

loglik=-nrow(X)*log(2*pi)/2-sum(w%*%delta)/2-sum(exp(-w%*%delta)*(Y-exp(X%*%beta))^2)/2

print(loglik)

w=cbind(rep(1,nrow(X)),exp(X%*%beta))

}








#p.126

#ARCH(多変量版) 交互勾配上昇法

data("AirPassengers")

Y=c(t(AirPassengers))

X=cbind(rep(1,length(Y)),c(1:length(Y)))

n=nrow(X)

beta=rep(1,ncol(X))

m=3

alpha=rep(1,m+1)

ite=10^5

eta=0.0001

for(l in 1:ite){

loglik=0;loglik_pre=loglik+1;

#alpha=rep(1,m+1)

while(sum(abs(loglik_pre-loglik))>10^(-2))({

loglik_pre=loglik  

u=Y-X%*%beta

dalpha=rep(0,length(alpha))

dbeta=rep(0,length(beta))

loglik=0

for(j in (m+1):(n)){

sub_u=u[(j-m):(j-1)]^2 

z=c(1,sub_u)

h=sum(alpha*c(1,sub_u))

dalpha_sub=z*(u[j]^2/h-1)/(2*h)

dalpha=dalpha+z*(u[j]^2/h-1)/(2*h)

loglik=loglik-log(ifelse(h>0,h,1))/2-u[j]^2/h

}

alpha=alpha+eta*dalpha

print(loglik)

})


loglik=0;loglik_pre=loglik+1;

while(sum(abs(loglik_pre-loglik))>10^(-2))({

loglik_pre=loglik  

u=Y-X%*%beta

dbeta=rep(0,length(beta))

loglik=0

for(j in (m+1):(n)){

sub_u=u[(j-m):(j-1)]^2 

z=c(1,sub_u)

h=sum(alpha*c(1,sub_u))

sub_val=-2*alpha[2:length(alpha)]*u[(j-m):(j-1)]

dbeta=dbeta+u[j]*X[j,]/h+apply(t(X[(j-m):(j-1),])*c(sub_val),1,sum)*(u[j]^2/h-1)/(2*h)

dbeta_sub=u[j]*X[j,]/h+apply(t(X[(j-m):(j-1),])*c(sub_val),1,sum)*(u[j]^2/h-1)/(2*h)

loglik=loglik-log(ifelse(h>0,h,1))/2-u[j]^2/h

}

beta=beta+eta*dbeta

#print(paste0(loglik,";",loglik_pre))

print(loglik)

})

}





#p.126

#ARCH(多変量版) BHHH法

#修正中

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56))

Y=data$y

X=cbind(rep(1,nrow(data)),as.matrix(data[,colnames(data) %in% c("x1","x2","x3")]))

n=nrow(X)

beta=rep(1,ncol(X))

m=3

alpha=rep(1,m+1)

ite=10^5

eta=0.0001

for(l in 1:ite){

loglik=0;loglik_pre=loglik+1;

#alpha=rep(1,m+1)

while(sum(abs(loglik_pre-loglik))>10^(-6))({

loglik_pre=loglik  
  
u=Y-X%*%beta

dalpha=rep(0,length(alpha))

dbeta=rep(0,length(beta))

ddalpha=array(0,dim=c(length(alpha),length(alpha)))

ddbeta=array(0,dim=c(length(beta),length(beta)))

loglik=0

for(j in (m+1):(n)){
  
sub_u=u[(j-m):(j-1)]^2 

z=c(1,sub_u)

h=sum(alpha*c(1,sub_u))

dalpha_sub=z*(u[j]^2/h-1)/(2*h)

dalpha=dalpha+z*(u[j]^2/h-1)/(2*h)

#sub_val=-2*alpha[2:length(alpha)]*u[(j-m):(j-1)]

#dbeta=dbeta+u[j]*X[j,]/h+apply(t(X[(j-m):(j-1),])*c(sub_val),1,sum)*(u[j]^2/h-1)/(2*h)

#dbeta_sub=u[j]*X[j,]/h+apply(t(X[(j-m):(j-1),])*c(sub_val),1,sum)*(u[j]^2/h-1)/(2*h)

ddalpha=ddalpha+t(t(c(dalpha_sub)))%*%t(c(dalpha_sub))

#ddbeta=ddbeta+t(t(c(dbeta_sub)))%*%t(c(dbeta_sub))

loglik=loglik-log(ifelse(h>0,h,1))/2-u[j]^2/h

}

alpha=alpha+eta*solve(ddalpha+diag(1,length(alpha)))%*%dalpha

#print(paste0(loglik,";",loglik_pre))

print(loglik)

})


loglik=0;loglik_pre=loglik+1;

#beta=rep(1,ncol(X))

while(sum(abs(loglik_pre-loglik))>10^(-6))({

loglik_pre=loglik  
  
u=Y-X%*%beta

#dalpha=rep(0,length(alpha))

dbeta=rep(0,length(beta))

#ddalpha=array(0,dim=c(length(alpha),length(alpha)))

ddbeta=array(0,dim=c(length(beta),length(beta)))

loglik=0

for(j in (m+1):(n)){
  
sub_u=u[(j-m):(j-1)]^2 

z=c(1,sub_u)

h=sum(alpha*c(1,sub_u))

#dalpha_sub=z*(u[j]^2/h-1)/(2*h)

#dalpha=dalpha+z*(u[j]^2/h-1)/(2*h)

sub_val=-2*alpha[2:length(alpha)]*u[(j-m):(j-1)]

dbeta=dbeta+u[j]*X[j,]/h+apply(t(X[(j-m):(j-1),])*c(sub_val),1,sum)*(u[j]^2/h-1)/(2*h)

dbeta_sub=u[j]*X[j,]/h+apply(t(X[(j-m):(j-1),])*c(sub_val),1,sum)*(u[j]^2/h-1)/(2*h)

#ddalpha=ddalpha+t(t(c(dalpha_sub)))%*%t(c(dalpha_sub))

ddbeta=ddbeta+t(t(c(dbeta_sub)))%*%t(c(dbeta_sub))

loglik=loglik-log(ifelse(h>0,h,1))/2-u[j]^2/h

}

beta=beta+eta*solve(ddbeta+diag(1,length(beta)))%*%dbeta

#print(paste0(loglik,";",loglik_pre))

print(loglik)

})

}






#p.164 Gauss法

data("AirPassengers")

Y=c(AirPassengers)

X=cumsum(sample(1:10,length(Y),replace = T))

XX=X[2:length(X)];XX_past=X[1:(length(X)-1)]

YY=Y[2:length(Y)];YY_past=Y[1:(length(Y)-1)]



f=function(x){

return(x[1]+x[2]*X^x[3])

}


h=0.001

ite=10^4

sita=rep(0.1,3)

eta=10^(-6)

for(l in 1:ite){

sita_pre=sita

Z=array(0,dim=c(length(X),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)%*%(Y-f(sita))

sita=c(sita+eta*beta)

print(sum(abs(sita_pre-sita)))

}








#ハウタッカー・テイラーモデル　p.240

#p.164 Gauss法

data("AirPassengers")

Y=c(AirPassengers)

X=cumsum(sample(1:10,length(Y),replace = T))

XX=X[2:length(X)];XX_past=X[1:(length(X)-1)]

YY=Y[2:length(Y)];YY_past=Y[1:(length(Y)-1)]



f=function(x){

ep=1-(x[2]-x[4])/2

return(x[1]*x[4]/ep+YY_past*(1+(x[2]-x[4])/2)/ep+x[3]*(1+x[4]/2)*(XX-XX_past)/ep+x[3]*x[4]*XX_past/ep)

}


h=0.001

ite=10^8

sita=rep(0.1,4)

eta=10^(-6)

for(l in 1:ite){

sita_pre=sita

Z=array(0,dim=c(length(XX),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)%*%(YY-f(sita))

sita=c(sita+eta*beta)

print(sum(abs(YY-f(sita))))

}






#p.241 (5)バーグストロムモデル

data("AirPassengers")

Y=c(AirPassengers)

X=cumsum(sample(1:10,length(Y),replace = T))

mat=c()

for(j in 2:length(Y)){
  
vec=c(1,Y[j-1],X[j]+X[j-1])
 
mat=c(mat,vec)
 
}

mat=t(matrix(mat,nrow=3))

y=Y[2:length(Y)]

A=solve(t(mat)%*%mat)%*%t(mat)%*%y

sita=2*(1-A[2])/(1+A[2])

alpha=A[1]/(1-A[2])

beta=2*A[3]/(1-A[2])

#A0

sita*alpha/(1+sita/2)

#A1

(1-sita/2)/(1+sita/2)

#A2

(sita*beta/2)/(1+sita/2)

#A3

sita*beta/(1+sita/2)

cbind(Y,mat%*%A)



#オリジナル　作成日：20210610

#自己回帰過程における予測子exp(X%*%beta)の場合のラッソ回帰

#勾配上昇法

data("AirPassengers")

p=3

data=data.frame(no=1:length(AirPassengers),values=c(AirPassengers))

r=10;dim=1

values=log(data$values)

X=array(0,dim=c((nrow(data)-r),dim*r))

Y=c()

for(j in 1:(nrow(data)-r)){

vec=c()

for(i in 1:dim){

vec=c(vec,values[j:(j+r-1)]^i)  

}

X[j,]=vec

Y=c(Y,(values[r+j]))

}



beta=rep(0.01,ncol(X))

eta=1/eigen(t(X)%*%c(exp(X%*%beta))%*%t(t(X)%*%c(exp(X%*%beta))))$values[1]

ite=10^7

ep=10^(-10)

for(l in 1:ite){

beta=beta+eta*t(X)%*%(exp(X%*%beta)*(Y-exp(X%*%beta)))

beta=ifelse(beta>ep,beta-ep,ifelse(beta<-ep,beta+ep,0))

cost=sum((Y-exp(X%*%beta))^2)+ep*sum(abs(beta))

print(cost)

#print(beta)

}
2
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
2
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?