xm = cbind(xm,prod((i-der+1):i)*x^(i-der))
}
else{
xm = cbind(xm,x^(i-der))
}
}
}
if(nknots > 0){
for(i in 1:(nknots)){
if((degree-der+1)<=degree){
xm = cbind(xm,prod((degree-der+1):degree) * (x-knots[i])^(degree-der)*(x > knots[i])) ;
}
else{
xm = cbind(xm,(x-knots[i])^(degree-der)*(x > knots[i])) ;
}
}
}
return(xm)
}
#简单调试
# n=600;
# t=runif(n,0,1)
# nknots=5;
# degree=4;
# knots = quantileknots(t,nknots);
# xm=basis_fun(t,degree,knots);
# Y=sin(2*pi*t)+rnorm(n,0,1);
#
# Beta=solve(t(xm)%*%xm)%*%t(xm)%*%Y;
# plot(t,sin(2*pi*t))
# lines(t,xm%*%Beta,col="red",type="p")
num=55
data=read.csv("wuhan.csv",header = T)
Yi=diff(data[,1])
x=as.matrix(data[-1,2:6])
x[,1]=x[,1]/10000
x[,3]=x[,3]/10000
Yi=Yi[1:num]
x=x[1:num,]
n=length(Yi)
Ui=seq(1,num)
nknots=c(2,4,3)
degree=c(2,5,4)
#对第一个变系数取样条基
nknots1=nknots[1]
degree1=degree[1]
knots1 = quantileknots(Ui,nknots1)
xm1=basis_fun(Ui,degree1,knots1)
x1=xm1*x[,1]
L1=ncol(x1)
#对第二个变系数取样条基
nknots2=nknots[2]
degree2=degree[2]
knots2 = quantileknots(Ui,nknots2)
xm2=basis_fun(Ui,degree2,knots2)
x2=xm2*x[,2]
L2=ncol(x2)
#对第三个变系数取样条基
nknots3=nknots[3]
degree3=degree[3]
knots3 = quantileknots(Ui,nknots3)
xm3=basis_fun(Ui,degree3,knots3)
x3=xm3*x[,3]
L3=ncol(x3)
#加权最小二乘求解系数
W=cbind(x1,x2,x3,x[,4],x[,5])
summ=glm(Yi~W,family=quasi(variance="mu",link="log"))
hatgam1=summ[[1]][2:(1+L1)]
hatgam2=summ[[1]][(1+L1+1):(1+L1+L2)]
hatgam3=summ[[1]][(1+L1+L2+1):(1+L1+L2+L3)]
hatbeta1Ui=xm1%*%hatgam1
hatbeta2Ui=xm2%*%hatgam2
hatbeta3Ui=xm3%*%hatgam3
hatbeta4=summ[[1]][1+L1+L2+L3+1]
hatbeta5=summ[[1]][1+L1+L2+L3+2]
#基于已有的系数，估计均值
lambda=numeric(n)
for(i in 1:n){
lambda[i]=exp(summ[[1]][1]+x[i,1]*hatbeta1Ui[i]+x[i,2]*hatbeta2Ui[i]+x[i,3]*hatbeta3Ui[i]+x[i,4]*hatbeta4+x[i,5]*hatbeta5)
}
write.csv(lambda,'lamhat.csv')
mean((Yi-lambda)^2)
par(mfrow=c(2,2))
plot(Ui,Yi,'l',col=1,xlim=c(min(Ui),max(Ui)),ylim=c(min(Yi),max(Yi)))
par(new=TRUE)
plot(Ui,lambda,type="l",pch=2,col=2,xlim=c(min(Ui),max(Ui)),ylim=c(min(Yi),max(Yi)),lwd=1,ann=F,xaxt='n',xlab='',ylab='')
plot(Ui,hatbeta1Ui,type="l",col=1,xlab="u",ylab="beta1",lwd=2)
plot(Ui,hatbeta2Ui,type="l",col=1,xlab="u",ylab="beta2",lwd=2)
plot(Ui,hatbeta3Ui,type="l",col=1,xlab="u",ylab="beta3",lwd=2)
#-------------计算m---------------
M=seq(1,30)
A=seq(30,5000,10)
MSE=matrix(0,length(A),length(M))
mstar=NULL
Yt=NULL
Lambda=lambda
for (l in (1:n)){
Yt[l]=data$Y[1]+sum(Lambda[1:l])
}
for (j in (1:length(A))){
a=A[j]
for (i in (1:length(M))){
Wk=NULL
Wk0=Lambda[1]+a
m=M[i]
print(m)
for (t in (1:(n-m))){
Wk[t]=Wk0+sum(Lambda[(m+1):(m+t)])
# Wk[t]=sum(Lambda[(m+1):(m+t)])+Wk0+sum(Lambda[1:t])
}
Wk=c(Wk0,Wk)
MSE[j,i]=mean((Yt[m:n]-Wk)^2)
}
mstar[j]=which(MSE[j,]==min(MSE[j,]))
}
which(MSE==min(MSE), arr.ind = TRUE)
#A=420,m=8
#——————————————————————Gamma--------------
M=seq(1,50)
a=420
m=8
Lambda=lambda
Wk=NULL
Wk0=Lambda[1]+a
for (t in (1:(n-m))){
Wk[t]=Wk0+sum(Lambda[(m+1):(m+t)])
# Wk[t]=sum(Lambda[(m+1):(m+t)])+Wk0+sum(Lambda[1:t])
}
Wk=c(Wk0,Wk)
Gamma=NULL
T=seq(m+1,n-m)
for (l in (1:length(T))){
tt=T[l]
Gamma[l]=(Wk[tt+1]-Wk[tt])/(Wk[tt]-Wk[tt-m])
}
gn=length(Gamma)
Xt=seq(1,gn)
Gamhat=NULL
h=1.3
for (i in (1:gn)){
Gamhat[i]=sum(Gamma*exp(-(Xt-Xt[i])^2/(2*h^2)))/sum(exp(-(Xt-Xt[i])^2/(2*h^2)))
}
#plot(Xt,Gamma,'l',xlim=c(0,gn),ylim=c(0,max(Gamma)))
#par(new=TRUE)
#plot(Xt,Gamhat,'l',xlim=c(0,gn),ylim=c(0,max(Gamma)),
#ann=F,col='red')
par(mfrow=c(1,2))
temp <- as.character(seq.Date(from = as.Date("2020/01/16"), by = "day", length.out = 56))
plot(T,Gamhat,'l',lty=1,xaxt='n',xlab='',ylab='Gamma')
axis(1,at=T,label=temp[T],las=3,cex.axis = 0.8,font.axis=10)
R=m*Gamhat
yt=expression(R[t])
plot(T,R,'l',lty=1,xaxt='n',xlab='',ylab=yt)
axis(1,at=T,label=temp[T],las=3,cex.axis = 0.8,font.axis=10)
write.csv(hatbeta1Ui,'beta1.csv')
write.csv(hatbeta2Ui,'beta2.csv')
write.csv(hatbeta3Ui,'beta3.csv')
write.csv(Gamhat,'Gamhat.csv')
rm(list=ls())
source('quantileknots.R')
source('basis_fun.R')
alldata=read.csv('wuhan.csv')
num=56
n=num-1
Times=1000
Ui=seq(1,n)
Bstbeta0=NULL
Bstbeta1=NULL
Bstbeta2=NULL
Bstbeta3=NULL
Bstbeta4=NULL
Bstbeta5=NULL
GAMMA=NULL
Yi=NULL
nknots1=2;degree1=2
nknots2=4;degree2=5
nknots3=3;degree3=4
data=alldata[(1:num),]
x=as.matrix(data[-1,2:6])
x[,1]=x[,1]/10000
x[,3]=x[,3]/10000
knots1 = quantileknots(Ui,nknots1)
xm1=basis_fun(Ui,degree1,knots1)
x1=xm1*x[,1]
L1=ncol(x1)
knots2 = quantileknots(Ui,nknots2)
xm2=basis_fun(Ui,degree2,knots2)
x2=xm2*x[,2]
L2=ncol(x2)
knots3 = quantileknots(Ui,nknots3)
xm3=basis_fun(Ui,degree3,knots3)
x3=xm3*x[,3]
L3=ncol(x3)
W=cbind(x1,x2,x3,x[,4],x[,5])
for (bst in (1:Times)){
print(bst)
lambda=as.matrix(read.csv('lamhat.csv'))[,2]
for (i in (1:n)){
Yi[i]=rpois(1,lambda[i])
}
summ=glm(Yi~W,family=quasi(variance="mu",link="log"))
hatgam1=summ[[1]][2:(1+L1)]
hatgam2=summ[[1]][(1+L1+1):(1+L1+L2)]
hatgam3=summ[[1]][(1+L1+L2+1):(1+L1+L2+L3)]
hatbeta1Ui=xm1%*%hatgam1
hatbeta2Ui=xm2%*%hatgam2
hatbeta3Ui=xm3%*%hatgam3
hatbeta4=summ[[1]][1+L1+L2+L3+1]
hatbeta5=summ[[1]][1+L1+L2+L3+2]
Bstbeta1=cbind(Bstbeta1,hatbeta1Ui)
Bstbeta2=cbind(Bstbeta2,hatbeta2Ui)
Bstbeta3=cbind(Bstbeta3,hatbeta3Ui)
Bstbeta0[bst]=summ[[1]][1]
Bstbeta4[bst]=hatbeta4
Bstbeta5[bst]=hatbeta5
lambda=numeric(n)
for(i in 1:n){
lambda[i]=exp(summ[[1]][1]+x[i,1]*hatbeta1Ui[i]+x[i,2]*hatbeta2Ui[i]+x[i,3]*hatbeta3Ui[i]+x[i,4]*hatbeta4+x[i,5]*hatbeta5)
}
a=420
m=8
Lambda=lambda
Wk=NULL
Wk0=Lambda[1]+a
for (t in (1:(n-m))){
Wk[t]=Wk0+sum(Lambda[(m+1):(m+t)])
# Wk[t]=sum(Lambda[(m+1):(m+t)])+Wk0+sum(Lambda[1:t])
}
Wk=c(Wk0,Wk)
Gamma=NULL
T=seq(m+1,n-m)
for (l in (1:length(T))){
tt=T[l]
Gamma[l]=(Wk[tt+1]-Wk[tt])/(Wk[tt]-Wk[tt-m])
}
GAMMA=cbind(GAMMA,Gamma)
}
GAMMA=rbind(matrix(0,m,1000),GAMMA,matrix(0,m,1000))
par(mfrow=c(1,2))
Ty=data$Y[-1]
Yhat=data$Y[1]+cumsum(Lambda)
plot(Ui,diff(data$Y),'l',col=1,xlim=c(min(Ui),max(Ui)),ylim=c(min(diff(data$Y)),max(diff(data$Y))),lwd=1,xlab='t',ylab='')
par(new=TRUE)
plot(Ui,lambda,type="l",lty=2,col=1,xlim=c(min(Ui),max(Ui)),ylim=c(min(diff(data$Y)),max(diff(data$Y))),lwd=1,ann=F,xaxt='n',yaxt='n')
yt=expression(lambda[paste('1t')])
title(main='(a)',font.main= 1, ylab=yt,ps=1,cex.main=1)
plot(Ui,Ty,'l',col=1,xlim=c(min(Ui),max(Ui)),ylim=c(min(Ty),max(Ty)),lwd=1,xlab='t',ylab='')
par(new=TRUE)
plot(Ui,Yhat,type="l",lty=2,col=1,xlim=c(min(Ui),max(Ui)),ylim=c(min(Ty),max(Ty)),lwd=1,ann=F,xaxt='n',yaxt='n')
yt=expression(Y[paste('1t')])
title(main='(b)',font.main= 1, ylab=yt,ps=1,cex.main=1)
beta1=as.matrix(read.csv('beta1.csv'))
beta2=as.matrix(read.csv('beta2.csv'))
beta3=as.matrix(read.csv('beta3.csv'))
b1l=NULL;b1u=NULL
b2l=NULL;b2u=NULL
b3l=NULL;b3u=NULL
gal=NULL;gau=NULL
for ( i in (1:n)){
s1=sort(Bstbeta1[i,])
b1l[i]=s1[1000*0.025];b1u[i]=s1[1000*0.975]
s2=sort(Bstbeta2[i,])
b2l[i]=s2[1000*0.025];b2u[i]=s2[1000*0.975]
s3=sort(Bstbeta3[i,])
b3l[i]=s3[1000*0.025];b3u[i]=s3[1000*0.975]
s4=sort(GAMMA[i,])
gal[i]=s4[1000*0.025];gau[i]=s4[1000*0.975]
}
gal=gal[9:47];gau=gau[9:47]
gn=length(gal)
Xt=seq(1,gn)
gau_hat=NULL;gal_hat=NULL
h=1.3
for (i in (1:gn)){
gau_hat[i]=sum(gau*exp(-(Xt-Xt[i])^2/(2*h^2)))/sum(exp(-(Xt-Xt[i])^2/(2*h^2)))
gal_hat[i]=sum(gal*exp(-(Xt-Xt[i])^2/(2*h^2)))/sum(exp(-(Xt-Xt[i])^2/(2*h^2)))
}
par(mfrow=c(2,3))
plot(Ui,beta1,type="l",ylim=c(min(b1l),max(b1u)),col=1,xlab="t",ylab=expression(beta[NN](t)),lwd=1)
par(new=TRUE)
plot(Ui,b1l,type="l",ylim=c(min(b1l),max(b1u)),lty=2,col=1,xaxt='n',yaxt='n',xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(Ui,b1u,type="l",ylim=c(min(b1l),max(b1u)),lty=2,col=1,xaxt='n',yaxt='n',xlab="",ylab="",lwd=1)
title(main='(a)',font.main= 1, ps=1,cex.main=1)
plot(Ui,beta2,type="l",ylim=c(min(b2l),max(b2u)),col=1,xlab="t",ylab=expression(beta[TN](t)),lwd=1)
par(new=TRUE)
plot(Ui,b2l,type="l",ylim=c(min(b2l),max(b2u)),lty=2,col=1,xaxt='n',yaxt='n',xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(Ui,b2u,type="l",ylim=c(min(b2l),max(b2u)),lty=2,col=1,xaxt='n',yaxt='n',xlab="",ylab="",lwd=1)
title(main='(b)',font.main= 1, ps=1,cex.main=1)
plot(Ui,beta3,type="l",ylim=c(min(b3l),max(b3u)),col=1,xlab="t",ylab=expression(beta[ND](t)),lwd=1)
par(new=TRUE)
plot(Ui,b3l,type="l",ylim=c(min(b3l),max(b3u)),lty=2,col=1,xaxt='n',yaxt='n',xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(Ui,b3u,type="l",ylim=c(min(b3l),max(b3u)),lty=2,col=1,xaxt='n',yaxt='n',xlab="",ylab="",lwd=1)
title(main='(c)',font.main= 1, ps=1,cex.main=1)
plot(Ui[20:n],beta1[20:n],type="l",ylim=c(min(b1l[20:n]),max(b1u[20:n])),col=1,xlab="t",ylab=expression(beta[NN](t)),lwd=1)
par(new=TRUE)
plot(Ui[20:n],b1l[20:n],type="l",ylim=c(min(b1l[20:n]),max(b1u[20:n])),lty=2,col=1,xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(Ui[20:n],b1u[20:n],type="l",ylim=c(min(b1l[20:n]),max(b1u[20:n])),lty=2,col=1,xlab="",ylab="",lwd=1)
title(main='(d)',font.main= 1, ps=1,cex.main=1)
plot(Ui[0:35],beta2[0:35],type="l",ylim=c(min(b2l[0:35]),max(b2u[0:35])),col=1,xlab="t",ylab=expression(beta[TN](t)),lwd=1)
par(new=TRUE)
plot(Ui[0:35],b2l[0:35],type="l",ylim=c(min(b2l[0:35]),max(b2u[0:35])),lty=2,col=1,xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(Ui[0:35],b2u[0:35],type="l",ylim=c(min(b2l[0:35]),max(b2u[0:35])),lty=2,col=1,xlab="",ylab="",lwd=1)
title(main='(e)',font.main= 1, ps=1,cex.main=1)
plot(Ui[15:35],beta3[15:35],type="l",ylim=c(min(b3l[15:35]),max(b3u[15:35])),col=1,xlab="t",ylab=expression(beta[ND](t)),lwd=1)
par(new=TRUE)
plot(Ui[15:35],b3l[15:35],type="l",ylim=c(min(b3l[15:35]),max(b3u[15:35])),lty=2,col=1,xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(Ui[15:35],b3u[15:35],type="l",ylim=c(min(b3l[15:35]),max(b3u[15:35])),lty=2,col=1,xlab="",ylab="",lwd=1)
title(main='(f)',font.main= 1, ps=1,cex.main=1)
# par(mfrow=c(1,2))
# plot(Ui[1:20],beta1[1:20],type="l",ylim=c(min(b1l[1:20]),max(b1u[1:20])),col=1,xlab="u",ylab="beta1",lwd=1)
# par(new=TRUE)
# plot(Ui[1:20],b1l[1:20],type="l",ylim=c(min(b1l[1:20]),max(b1u[1:20])),lty=3,col=2,xlab="",ylab="",lwd=1)
# par(new=TRUE)
# plot(Ui[1:20],b1u[1:20],type="l",ylim=c(min(b1l[1:20]),max(b1u[1:20])),lty=3,col=2,xlab="",ylab="",lwd=1)
#
# plot(Ui[20:n],beta1[20:n],type="l",ylim=c(min(b1l[20:n]),max(b1u[20:n])),col=1,xlab="u",ylab="beta1",lwd=1)
# par(new=TRUE)
# plot(Ui[20:n],b1l[20:n],type="l",ylim=c(min(b1l[20:n]),max(b1u[20:n])),lty=3,col=2,xlab="",ylab="",lwd=1)
# par(new=TRUE)
# plot(Ui[20:n],b1u[20:n],type="l",ylim=c(min(b1l[20:n]),max(b1u[20:n])),lty=3,col=2,xlab="",ylab="",lwd=1)
#
# par(mfrow=c(1,2))
# plot(Ui[1:35],beta2[1:35],type="l",ylim=c(min(b2l[1:35]),max(b2u[1:35])),col=1,xlab="u",ylab="beta2",lwd=1)
# par(new=TRUE)
# plot(Ui[1:35],b2l[1:35],type="l",ylim=c(min(b2l[1:35]),max(b2u[1:35])),lty=3,col=2,xlab="",ylab="",lwd=1)
# par(new=TRUE)
# plot(Ui[1:35],b2u[1:35],type="l",ylim=c(min(b2l[1:35]),max(b2u[1:35])),lty=3,col=2,xlab="",ylab="",lwd=1)
#
# plot(Ui[35:n],beta2[35:n],type="l",ylim=c(min(b2l[35:n]),max(b2u[35:n])),col=1,xlab="u",ylab="beta2",lwd=1)
# par(new=TRUE)
# plot(Ui[35:n],b2l[35:n],type="l",ylim=c(min(b2l[35:n]),max(b2u[35:n])),lty=3,col=2,xlab="",ylab="",lwd=1)
# par(new=TRUE)
# plot(Ui[35:n],b2u[35:n],type="l",ylim=c(min(b2l[35:n]),max(b2u[35:n])),lty=3,col=2,xlab="",ylab="",lwd=1)
#
# par(mfrow=c(1,3))
# plot(Ui[1:15],beta3[1:15],type="l",ylim=c(min(b3l[1:15]),max(b3u[1:15])),col=1,xlab="u",ylab="beta3",lwd=1)
# par(new=TRUE)
# plot(Ui[1:15],b3l[1:15],type="l",ylim=c(min(b3l[1:15]),max(b3u[1:15])),lty=3,col=2,xlab="",ylab="",lwd=1)
# par(new=TRUE)
# plot(Ui[1:15],b3u[1:15],type="l",ylim=c(min(b3l[1:15]),max(b3u[1:15])),lty=3,col=2,xlab="",ylab="",lwd=1)
#
# plot(Ui[15:35],beta3[15:35],type="l",ylim=c(min(b3l[15:35]),max(b3u[15:35])),col=1,xlab="u",ylab="beta3",lwd=1)
# par(new=TRUE)
# plot(Ui[15:35],b3l[15:35],type="l",ylim=c(min(b3l[15:35]),max(b3u[15:35])),lty=3,col=2,xlab="",ylab="",lwd=1)
# par(new=TRUE)
# plot(Ui[15:35],b3u[15:35],type="l",ylim=c(min(b3l[15:35]),max(b3u[15:35])),lty=3,col=2,xlab="",ylab="",lwd=1)
#
# plot(Ui[35:n],beta3[35:n],type="l",ylim=c(min(b3l[40:n]),max(b3u[40:n])),col=1,xlab="u",ylab="beta3",lwd=1)
# par(new=TRUE)
# plot(Ui[35:n],b3l[35:n],type="l",ylim=c(min(b3l[40:n]),max(b3u[40:n])),lty=3,col=2,xlab="",ylab="",lwd=1)
# par(new=TRUE)
# plot(Ui[35:n],b3u[35:n],type="l",ylim=c(min(b3l[40:n]),max(b3u[40:n])),lty=3,col=2,xlab="",ylab="",lwd=1)
Gamhat=as.matrix(read.csv('Gamhat.csv'))
par(mfrow=c(1,2))
T=seq(9,47)
temp <- as.character(seq.Date(from = as.Date("2020/01/16"), by = "day", length.out = 56))
plot(T,Gamhat,type="l",ylim=c(min(gal_hat),max(gau_hat)),lty=1,col=1,xaxt='n',xlab="",ylab="Gamma_1t",lwd=1)
par(new=TRUE)
plot(T,gal_hat,type="l",ylim=c(min(gal_hat),max(gau_hat)),lty=3,col=2,xaxt='n',xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(T,gau_hat,type="l",ylim=c(min(gal_hat),max(gau_hat)),lty=3,col=2,xaxt='n',xlab="",ylab="",lwd=1)
axis(1,at=T,label=temp[T],las=3,cex.axis = 0.8,font.axis=10)
title(main='(a)',font.main= 1, ps=1,cex.main=1)
T=seq(7,44)
temp <- as.character(seq.Date(from = as.Date("2020/01/28"), by = "day", length.out = 56))
plot(T,gamhat,type="l",ylim=c(min(new_gamhat_low),max(new_gamhat_upper)),lty=1,col=1,xaxt='n',xlab="",ylab="Gamma_Kt",lwd=1)
par(new=TRUE)
plot(T,new_gamhat_low,type="l",ylim=c(min(new_gamhat_low),max(new_gamhat_upper)),lty=3,col=2,xaxt='n',xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(T,new_gamhat_upper,type="l",ylim=c(min(new_gamhat_low),max(new_gamhat_upper)),lty=3,col=2,xaxt='n',xlab="",ylab="",lwd=1)
axis(1,at=T,label=temp[T],las=3,cex.axis = 0.8,font.axis=10)
title(main='(b)',font.main= 1, ps=1,cex.main=1)
par(mfrow=c(1,2))
T=seq(9,47)
temp <- as.character(seq.Date(from = as.Date("2020/01/16"), by = "day", length.out = 56))
plot(T,8*Gamhat,type="l",ylim=c(min(8*gal_hat),max(8*gau_hat)),lty=1,col=1,xaxt='n',xlab="",ylab=expression(R[paste('1t')]),lwd=1)
par(new=TRUE)
plot(T,8*gal_hat,type="l",ylim=c(min(8*gal_hat),max(8*gau_hat)),lty=2,col=1,xaxt='n',yaxt='n',xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(T,8*gau_hat,type="l",ylim=c(min(8*gal_hat),max(8*gau_hat)),lty=2,col=1,xaxt='n',yaxt='n',xlab="",ylab="",lwd=1)
abline(v=20,col='red',lty=2)
abline(v=36,col='red',lty=2)
axis(1,at=T,label=temp[T],las=3,cex.axis = 0.8,font.axis=10)
title(main='(a)',font.main= 1, ps=1,cex.main=1)
T=seq(7,44)
yt=expression(bar(R)[t])
temp <- as.character(seq.Date(from = as.Date("2020/01/28"), by = "day", length.out = 56))
plot(T,6*gamhat,type="l",ylim=c(min(6*new_gamhat_low),max(6*new_gamhat_upper)),lty=1,col=1,xaxt='n',xlab="",ylab=yt,lwd=1)
par(new=TRUE)
plot(T,6*new_gamhat_low,type="l",ylim=c(min(6*new_gamhat_low),max(6*new_gamhat_upper)),lty=2,col=1,xaxt='n',yaxt='n',xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(T,6*new_gamhat_upper,type="l",ylim=c(min(6*new_gamhat_low),max(6*new_gamhat_upper)),lty=2,col=1,xaxt='n',yaxt='n',xlab="",ylab="",lwd=1)
abline(v=9,col='red',lty=2)
abline(v=21,col='red',lty=2)
axis(1,at=T,label=temp[T],las=3,cex.axis = 0.8,font.axis=10)
title(main='(b)',font.main= 1, ps=1,cex.main=1)
par(mfrow=c(1,2))
Ty=data$Y[-1]
Yhat=data$Y[1]+cumsum(Lambda)
plot(Ui,diff(data$Y),'l',col=1,xlim=c(min(Ui),max(Ui)),ylim=c(min(diff(data$Y)),max(diff(data$Y))),lwd=1,xlab='t',ylab='')
par(new=TRUE)
plot(Ui,lambda,type="l",lty=2,col=1,xlim=c(min(Ui),max(Ui)),ylim=c(min(diff(data$Y)),max(diff(data$Y))),lwd=1,ann=F,xaxt='n',yaxt='n')
yt=expression(lambda[paste('1t')])
title(main='(a)',font.main= 1, ylab=yt,ps=1,cex.main=1)
plot(Ui,Ty,'l',col=1,xlim=c(min(Ui),max(Ui)),ylim=c(min(Ty),max(Ty)),lwd=1,xlab='t',ylab='')
par(new=TRUE)
plot(Ui,Yhat,type="l",lty=2,col=1,xlim=c(min(Ui),max(Ui)),ylim=c(min(Ty),max(Ty)),lwd=1,ann=F,xaxt='n',yaxt='n')
yt=expression(Y[paste('1t')])
title(main='(b)',font.main= 1, ylab=yt,ps=1,cex.main=1)
beta1=as.matrix(read.csv('beta1.csv'))
par(mfrow=c(2,3))
plot(Ui,beta1,type="l",ylim=c(min(b1l),max(b1u)),col=1,xlab="t",ylab=expression(beta[NN](t)),lwd=1)
beta1=as.matrix(read.csv('beta1.csv'))[,2]
beta2=as.matrix(read.csv('beta2.csv'))[,2]
beta3=as.matrix(read.csv('beta3.csv'))[,2]
plot(Ui,beta1,type="l",ylim=c(min(b1l),max(b1u)),col=1,xlab="t",ylab=expression(beta[NN](t)),lwd=1)
par(mfrow=c(2,3))
plot(Ui,beta1,type="l",ylim=c(min(b1l),max(b1u)),col=1,xlab="t",ylab=expression(beta[NN](t)),lwd=1)
par(new=TRUE)
plot(Ui,b1l,type="l",ylim=c(min(b1l),max(b1u)),lty=2,col=1,xaxt='n',yaxt='n',xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(Ui,b1u,type="l",ylim=c(min(b1l),max(b1u)),lty=2,col=1,xaxt='n',yaxt='n',xlab="",ylab="",lwd=1)
title(main='(a)',font.main= 1, ps=1,cex.main=1)
plot(Ui,beta2,type="l",ylim=c(min(b2l),max(b2u)),col=1,xlab="t",ylab=expression(beta[TN](t)),lwd=1)
par(new=TRUE)
plot(Ui,b2l,type="l",ylim=c(min(b2l),max(b2u)),lty=2,col=1,xaxt='n',yaxt='n',xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(Ui,b2u,type="l",ylim=c(min(b2l),max(b2u)),lty=2,col=1,xaxt='n',yaxt='n',xlab="",ylab="",lwd=1)
title(main='(b)',font.main= 1, ps=1,cex.main=1)
plot(Ui,beta3,type="l",ylim=c(min(b3l),max(b3u)),col=1,xlab="t",ylab=expression(beta[ND](t)),lwd=1)
par(new=TRUE)
plot(Ui,b3l,type="l",ylim=c(min(b3l),max(b3u)),lty=2,col=1,xaxt='n',yaxt='n',xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(Ui,b3u,type="l",ylim=c(min(b3l),max(b3u)),lty=2,col=1,xaxt='n',yaxt='n',xlab="",ylab="",lwd=1)
title(main='(c)',font.main= 1, ps=1,cex.main=1)
plot(Ui[20:n],beta1[20:n],type="l",ylim=c(min(b1l[20:n]),max(b1u[20:n])),col=1,xlab="t",ylab=expression(beta[NN](t)),lwd=1)
par(new=TRUE)
plot(Ui[20:n],b1l[20:n],type="l",ylim=c(min(b1l[20:n]),max(b1u[20:n])),lty=2,col=1,xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(Ui[20:n],b1u[20:n],type="l",ylim=c(min(b1l[20:n]),max(b1u[20:n])),lty=2,col=1,xlab="",ylab="",lwd=1)
title(main='(d)',font.main= 1, ps=1,cex.main=1)
plot(Ui[0:35],beta2[0:35],type="l",ylim=c(min(b2l[0:35]),max(b2u[0:35])),col=1,xlab="t",ylab=expression(beta[TN](t)),lwd=1)
par(new=TRUE)
plot(Ui[0:35],b2l[0:35],type="l",ylim=c(min(b2l[0:35]),max(b2u[0:35])),lty=2,col=1,xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(Ui[0:35],b2u[0:35],type="l",ylim=c(min(b2l[0:35]),max(b2u[0:35])),lty=2,col=1,xlab="",ylab="",lwd=1)
title(main='(e)',font.main= 1, ps=1,cex.main=1)
plot(Ui[15:35],beta3[15:35],type="l",ylim=c(min(b3l[15:35]),max(b3u[15:35])),col=1,xlab="t",ylab=expression(beta[ND](t)),lwd=1)
par(new=TRUE)
plot(Ui[15:35],b3l[15:35],type="l",ylim=c(min(b3l[15:35]),max(b3u[15:35])),lty=2,col=1,xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(Ui[15:35],b3u[15:35],type="l",ylim=c(min(b3l[15:35]),max(b3u[15:35])),lty=2,col=1,xlab="",ylab="",lwd=1)
title(main='(f)',font.main= 1, ps=1,cex.main=1)
Gamhat=as.matrix(read.csv('Gamhat.csv'))
par(mfrow=c(1,2))
T=seq(9,47)
temp <- as.character(seq.Date(from = as.Date("2020/01/16"), by = "day", length.out = 56))
plot(T,Gamhat,type="l",ylim=c(min(gal_hat),max(gau_hat)),lty=1,col=1,xaxt='n',xlab="",ylab="Gamma_1t",lwd=1)
par(new=TRUE)
plot(T,gal_hat,type="l",ylim=c(min(gal_hat),max(gau_hat)),lty=3,col=2,xaxt='n',xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(T,gau_hat,type="l",ylim=c(min(gal_hat),max(gau_hat)),lty=3,col=2,xaxt='n',xlab="",ylab="",lwd=1)
axis(1,at=T,label=temp[T],las=3,cex.axis = 0.8,font.axis=10)
title(main='(a)',font.main= 1, ps=1,cex.main=1)
T=seq(7,44)
temp <- as.character(seq.Date(from = as.Date("2020/01/28"), by = "day", length.out = 56))
plot(T,gamhat,type="l",ylim=c(min(new_gamhat_low),max(new_gamhat_upper)),lty=1,col=1,xaxt='n',xlab="",ylab="Gamma_Kt",lwd=1)
par(new=TRUE)
plot(T,new_gamhat_low,type="l",ylim=c(min(new_gamhat_low),max(new_gamhat_upper)),lty=3,col=2,xaxt='n',xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(T,new_gamhat_upper,type="l",ylim=c(min(new_gamhat_low),max(new_gamhat_upper)),lty=3,col=2,xaxt='n',xlab="",ylab="",lwd=1)
axis(1,at=T,label=temp[T],las=3,cex.axis = 0.8,font.axis=10)
title(main='(b)',font.main= 1, ps=1,cex.main=1)
Gamhat=as.matrix(read.csv('Gamhat.csv'))[,2]
par(mfrow=c(1,2))
T=seq(9,47)
temp <- as.character(seq.Date(from = as.Date("2020/01/16"), by = "day", length.out = 56))
plot(T,Gamhat,type="l",ylim=c(min(gal_hat),max(gau_hat)),lty=1,col=1,xaxt='n',xlab="",ylab="Gamma_1t",lwd=1)
par(new=TRUE)
plot(T,gal_hat,type="l",ylim=c(min(gal_hat),max(gau_hat)),lty=3,col=2,xaxt='n',xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(T,gau_hat,type="l",ylim=c(min(gal_hat),max(gau_hat)),lty=3,col=2,xaxt='n',xlab="",ylab="",lwd=1)
axis(1,at=T,label=temp[T],las=3,cex.axis = 0.8,font.axis=10)
title(main='(a)',font.main= 1, ps=1,cex.main=1)
T=seq(7,44)
temp <- as.character(seq.Date(from = as.Date("2020/01/28"), by = "day", length.out = 56))
plot(T,gamhat,type="l",ylim=c(min(new_gamhat_low),max(new_gamhat_upper)),lty=1,col=1,xaxt='n',xlab="",ylab="Gamma_Kt",lwd=1)
par(new=TRUE)
plot(T,new_gamhat_low,type="l",ylim=c(min(new_gamhat_low),max(new_gamhat_upper)),lty=3,col=2,xaxt='n',xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(T,new_gamhat_upper,type="l",ylim=c(min(new_gamhat_low),max(new_gamhat_upper)),lty=3,col=2,xaxt='n',xlab="",ylab="",lwd=1)
axis(1,at=T,label=temp[T],las=3,cex.axis = 0.8,font.axis=10)
title(main='(b)',font.main= 1, ps=1,cex.main=1)
par(mfrow=c(1,2))
T=seq(9,47)
temp <- as.character(seq.Date(from = as.Date("2020/01/16"), by = "day", length.out = 56))
plot(T,Gamhat,type="l",ylim=c(min(gal_hat),max(gau_hat)),lty=1,col=1,xaxt='n',xlab="",ylab="Gamma_1t",lwd=1)
par(new=TRUE)
plot(T,gal_hat,type="l",ylim=c(min(gal_hat),max(gau_hat)),lty=3,col=2,xaxt='n',xlab="",ylab="",lwd=1)
par(new=TRUE)
plot(T,gau_hat,type="l",ylim=c(min(gal_hat),max(gau_hat)),lty=3,col=2,xaxt='n',xlab="",ylab="",lwd=1)
axis(1,at=T,label=temp[T],las=3,cex.axis = 0.8,font.axis=10)
title(main='(a)',font.main= 1, ps=1,cex.main=1)
T=seq(7,44)
temp <- as.character(seq.Date(from = as.Date("2020/01/28"), by = "day", length.out = 56))
plot(T,gamhat,type="l",ylim=c(min(new_gamhat_low),max(new_gamhat_upper)),lty=1,col=1,xaxt='n',xlab="",ylab="Gamma_Kt",lwd=1)
