set.seed(9)
n.1 <- 200
y.1 <- numeric(n.1)
sigma2.1 <- numeric(n.1)
sigma2.1[1] <- 0.05
z <- rnorm(1,sd=sqrt(sigma2[1]))
y.1[1]<-z^2


for(i in 2:n.1){
    sigma2.1[i] <- 0.2*z^2+0.75*sigma2.1[i-1]+0.05
    z <- rnorm(1,sd=sqrt(sigma2.1[i]))
    y.1[i]<-z^2
}

plot(1:n.1,y.1,pch=19,xlab="Time",ylab="Y(t)")
lines(1:n.1,sigma2.1,col=2,lwd=3)
abline(a=0.05,b=0,col=3,lwd=3)
abline(a=5,b=0,col=4,lwd=3)
legend("topleft",legend=c("Statiscian","Pessimist","Opsinist"),lty=1,lwd=3,col=2:4)

## Now with large n
n <- 200000
y <- numeric(n)
sigma2 <- numeric(n)
sigma2[1] <- 0.05
z <- rnorm(1,sd=sqrt(sigma2[1]))
y[1]<-z^2


for(i in 2:n){
    sigma2[i] <- 0.2*z^2+0.75*sigma2[i-1]+0.05
    z <- rnorm(1,sd=sqrt(sigma2[i]))
    y[i]<-z^2
}



## Scoreing rules
SE <- function(x,y){mean((x-y)^2)}
AE <- function(x,y){mean(abs(x-y))}
APE <- function(x,y)(mean(abs((x-y)/y)))
RE <- function(x,y)(mean(abs((x-y)/x)))

## 3 competitors
## Statistician: y_t = sigma2_t
## Optimist: y_t = 5
## Pessimist y_t = 0.05


## SE
se <- c("stat"=SE(sigma2,y),"optimist"=SE(5,y),"pessimist"=SE(0.05,y))
se


## AE
ae <- c("stat"=AE(sigma2,y),"optimist"=AE(5,y),"pessimist"=AE(0.05,y))
ae

## APE
ape <- c("stat"=APE(sigma2,y),"optimist"=APE(5,y),"pessimist"=APE(0.05,y))
ape/n


## re
re <- c("stat"=RE(sigma2,y),"optimist"=RE(5,y),"pessimist"=RE(0.05,y))
re


##################################################
## New competitor Mr. Bayes, Mr Bayes use knowledge of the predictive distribution
##################################################
yh.se <- sigma2
yh.ae <- qchisq(0.5,df=1) * sigma2
yh.ape <- 1e-10 ## i.e a very small number
yh.re <- qchisq(0.5,df=3) * sigma2

## SE
se <- c(se,"bayes" = SE(yh.se,y))
se


## AE
ae <- c(ae,"bayes" = AE(yh.ae,y))
ae


## APE
ape <- c(ape,"bayes" = APE(yh.ape,y))
ape


## RE
re <- c(re,"bayes" = RE(yh.re,y))
re


##################################################
## illustratinon MR. Bayes foresacts
yh.se <- sigma2.1
yh.ae <- qchisq(0.5,df=1) * sigma2.1
yh.ape <- 1e-10 ## i.e a very small number
yh.re <- qchisq(0.5,df=3) * sigma2.1


plot(1:n.1,y.1,pch=19,xlab="Time",ylab="Y(t)")
lines(1:n.1,yh.se,col=2,lwd=3)
lines(1:n.1,yh.ae,col=3,lwd=3)
abline(a=yh.ape,b=0,col=4,lwd=3)
lines(1:n.1,yh.re,col=5,lwd=3)
legend("topleft",legend=c("SE","AE","APE","RE"),lty=1,lwd=3,col=2:5)


##################################################
## PIT
x <- seq(0,30,length=300)
par(mfrow=c(2,2))
z.1 <- pchisq(y.1/sigma2.1,df=1)
z <- pchisq(y/sigma2,df=1)
hist(y.1/sigma2.1,prob=TRUE)
lines(x,dchisq(x,df=1),col=4,lwd=3)
hist(z.1)
hist(y/sigma2,prob=TRUE)
lines(x,dchisq(x,df=1),col=4,lwd=3)
hist(z)

##################################################
## end
##################################################