COVID19
COVID19 copied to clipboard
Kalman Filter R code reference
Seems like the Kalman filter has been implemented in R, is it possible to give reference to the mathematical formulae used while implementing it.
Especially in the following code
date[nrow(date) + 1,1] <-all[nrow(all),1]+1
pred_all<-NULL
for (n in 2:ncol(all)-1) {
Y<-ts(data = all[n+1], start = 1, end =nrow(all)+1)
sig_w<-0.01
w<-sig_w*randn(1,100) # acceleration which denotes the fluctuation (Q/R) rnorm(100, mean = 0, sd = 1)
sig_v<-0.01
v<-sig_v*randn(1,100)
t<-0.45
phi<-matrix(c(1,0,t,1),2,2)
gama<-matrix(c(0.5*t^2,t),2,1)
H<-matrix(c(1,0),1,2)
#Kalman
x0_0<-p0_0<-matrix(c(0,0),2,1)
p0_0<-matrix(c(1,0,0,1),2,2)
Q<-0.01
R<-0.01
X<-NULL
X2<-NULL
pred<-NULL
for (i in 0:nrow(all)) {
namp <-paste("p", i+1,"_",i, sep = "")
assign(namp, phi%*%(get(paste("p", i,"_",i, sep = "")))%*%t(phi)+gama%*%Q%*%t(gama))
namk <- paste("k", i+1, sep = "")
assign(namk,get(paste("p", i+1,"_",i, sep = ""))%*%t(H)%*%(1/(H%*%get(paste("p", i+1,"_",i, sep = ""))%*%t(H)+R)))
namx <- paste("x", i+1,"_",i, sep = "")
assign(namx,phi%*%get(paste("x", i,"_",i, sep = "")))
namE <- paste("E", i+1, sep = "")
assign(namE,Y[i+1]-H%*%get(paste("x", i+1,"_",i, sep = "")))
namx2 <- paste("x", i+1,"_",i+1, sep = "")
assign(namx2,get(paste("x", i+1,"_",i, sep = ""))+get(paste("k", i+1, sep = ""))%*%get(paste("E", i+1, sep = "")))
namp2 <- paste("p", i+1,"_",i+1, sep = "")
assign(namp2,(p0_0-get(paste("k", i+1, sep = ""))%*%H)%*%get(paste("p", i+1,"_",i, sep = "")))
X<-rbind(X,get(paste("x", i+1,"_",i,sep = ""))[1])
X2<-rbind(X2,get(paste("x", i+1,"_",i,sep = ""))[2])
if(i>2){
remove(list=(paste("p", i-1,"_",i-2, sep = "")))
remove(list=(paste("k", i-1, sep = "")))
remove(list=(paste("E", i-1, sep = "")))
remove(list=(paste("p", i-2,"_",i-2, sep = "")))
remove(list=(paste("x", i-1,"_",i-2, sep = "")))
remove(list=(paste("x", i-2,"_",i-2, sep = "")))}
}
pred<-NULL
pred<-cbind(Y,X,round(X2,4))
pred<-as.data.frame(pred)
pred$region<-colnames(all[,n+1])
pred$date<-date$date
pred$actual<-rbind(0,(cbind(pred[2:nrow(pred),1])/pred[1:nrow(pred)-1,1]-1)*100)
pred$predict<-rbind(0,(cbind(pred[2:nrow(pred),2])/pred[1:nrow(pred)-1,2]-1)*100)
pred$pred_rate<-(pred$X/pred$Y-1)*100
pred$X2_change<-rbind(0,(cbind(pred[2:nrow(pred),3]-pred[1:nrow(pred)-1,3])))
pred_all<-rbind(pred_all,pred)
}
Thank you