Theory and Practice of Recursive Identification, Ljung and soderstrom

#p.17 The least squares Method

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),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),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22))

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

Y=data$y*100

sita=rep(0,ncol(X))

a=sample(1:9,nrow(X),replace=T)

X=(log(X^2))

R=array(0,dim=c(ncol(X),ncol(X)))

vec=rep(0,ncol(X))

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

R=R+(a[j]*t(t(c(X[j,])))%*%t(c(X[j,])))

R_inv=svd(R)$u%*%diag(1/svd(R)$d)%*%t(svd(R)$v)

vec=vec+a[j]*X[j,]*Y[j]

sita=R_inv%*%vec  

}

W=diag(a)

solve(t(X)%*%W%*%X)%*%t(X)%*%W%*%Y

#p.20 The least squares Method(2.21)

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),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),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22))

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

Y=data$y*100

sita=rep(0,ncol(X))

a=sample(1:9,nrow(X),replace=T)

X=(log(X^2))

R=array(0,dim=c(ncol(X),ncol(X)))

vec=rep(0,ncol(X))

P=diag(rep(10000,ncol(X)))

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

R=R+(a[j]*t(t(c(X[j,])))%*%t(c(X[j,])))

vec=vec+solve(P)%*%sita+a[j]*X[j,]*Y[j]

sita=solve(solve(P)+R)%*%vec  

}

W=diag(a)

solve(t(X)%*%W%*%X)%*%t(X)%*%W%*%Y

#p.20 The recursive least squares Method(2.19a),(2.19b),(2.19c)

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),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),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22))

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

Y=data$y*100

sita=rep(0,ncol(X))

a=sample(1:9,nrow(X),replace=T)

X=(log(X^2))

R=array(0,dim=c(ncol(X),ncol(X)))

vec=rep(0,ncol(X))

P=diag(rep(100000,ncol(X)))

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

L=P%*%X[j,]/as.numeric(1/a[j] + t(c(X[j, ])) %*% P %*% t(t(c(X[j, ]))))

P=P-P%*%t(t(c(X[j,])))%*%t(c(X[j,]))%*%P/as.numeric(1/a[j] + t(c(X[j, ])) %*% P %*% t(t(c(X[j, ]))))

sita=sita+L%*%(Y[j]-sum(X[j,]*sita))

}

W=diag(a)

solve(t(X)%*%W%*%X)%*%t(X)%*%W%*%Y

# 2.2.3 A Recursive Prediction Error Method(ARMAX)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82))

Y=data$y;u=data$x1

a=0.1;b=0.1;c=0.1;e=0.1

param=c(a,b,c)

R=array(0,dim=c(3,3))

e_vec=c()

for(j in 1:length(Y)){

pthi=c(-Y[j],u[j],e)

R=R+t(t(c(pthi)))%*%t(c(pthi))

pre=sum(pthi*param)

if(j>1){

e=Y[j]+sum(param*c(Y[j-1],-u[j-1],-1))

}

e_vec=c(e_vec,e)

param=param+solve(R+diag(rep(0.0001,length(param))))%*%pthi*(Y[j]-pre)

}

N=nrow(data);e_vec=e_vec[-N]

predict=cbind(-Y[1:(N-1)],u[1:(N-1)],e_vec)%*%param

#p.87 3.3.5 Optimal Choice of Weights in Quadratic Criteria

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),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),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22))

X=as.matrix(data[,colnames(data) %in% c("x1","x2","x3")])/100

Y=data$y*100

n=nrow(data);r=1

param=rep(1,ncol(X))

lam=diag(rep(1,n))

R=array(0,dim=c(ncol(X),ncol(X)))

for(j in 1:n){

r=1/j

pre=X%*%param

epsilon=Y-pre

lam=lam+(t(t(c(epsilon)))%*%t(c(epsilon))-lam)*r

R=R+r*(t(X)%*%solve(lam+diag(rep(0.00001,nrow(X))))%*%X-R)

sita=sita-r*solve(R+diag(rep(0.00001,ncol(X))))%*%t(X)%*%solve(lam+diag(rep(0.00001,nrow(X))))%*%epsilon

cost=t(epsilon)%*%solve(lam+diag(rep(0.00001,nrow(X))))%*%epsilon

print(cost)

}