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

```{r}


#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)

```


```{r}



#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))

}





```

```{r}


#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)

})

}



```


```{r}

#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)

})

}



```

```{r}


#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)))

}





```


```{r}


#ハウタッカー・テイラーモデル　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))))

}




```

```{r}

#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)

```

```{r}

#オリジナル　作成日：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)

}





```