###
### Artificial PET data
###
x <- 1:128
f12 <- function(x,y){
2*((1.5*(x-64)^2+(y-64)^2<3025)) +
((x-64)^2+(y-104)^2<16)-((x-92)^2+(y-84)^2<25)+
((x-78)^2+(y-84)^2<25)-((x-64)^2+(y-84)^2<36)+
((x-50)^2+(y-84)^2<36)-((x-36)^2+(y-84)^2<25)+
((x-92)^2+(y-64)^2<25)-((x-78)^2+(y-64)^2<16)+
((x-64)^2+(y-64)^2<16)-((x-50)^2+(y-64)^2<25)+
((x-36)^2+(y-64)^2<25)-((x-92)^2+(y-44)^2<36)+
((x-78)^2+(y-44)^2<36)-((x-64)^2+(y-44)^2<25)+
((x-50)^2+(y-44)^2<25)-((x-36)^2+(y-44)^2<16)+
((x-64)^2+(y-24)^2<16)
}
u0 <- 2*outer(x,x,"f12")
set.seed(1)
y <- matrix(rpois(u0,u0),128,128)
# use larger hmax for good results
yhat <- laws(y,model="Poisson",hmax=4)$theta
par(mfrow=c(1,3),mar=c(3,3,3,.5),mgp=c(2,1,0))
image(y,col=gray((0:255)/255))
title("Observed image")
image(yhat,col=gray((0:255)/255))
title("AWS-Reconstruction")
image(u0,col=gray((0:255)/255))
title("True image")
rm(u0,y,yhat,x)
###
### VOLATITILTY ESTIMATION
###
### artificial example
###
sigma <- c(rep(1,125),rep(2,125),rep(.5,125),rep(1,125))
set.seed(1)
x <- rnorm(sigma,0,sigma)
tmpa <- laws(x,model="Volatility",u=sigma,graph=TRUE,hmax=250)
tmps <- laws(x,model="Volatility",u=sigma,hmax=250,symmetric=TRUE)
par(mfrow=c(1,1),mar=c(3,3,3,1))
plot(abs(x),col=3,xlab="time t",ylab=expression(abs( R[t] )))
lines(tmpa$theta,col=1,lwd=2)
lines(tmps$theta,col=1,lwd=2,lty=2)
lines(sigma,col=2,lwd=2,lty=3)
legend(350,5.5,c("asymmetric AWS","symmetric AWS","true sigma"),
lwd=c(2,2,2),lty=c(1,2,3),col=c(1,1,2))
title(expression(paste("Estimates of ",sigma(t))))
rm(tmpa,tmps,sigma,x)Run the code above in your browser using DataLab