## com303a of september 2026
#  re stk4011 Ch3 and mse functions
#  simple model: Y is N(\theta,1)
#  three estimators: Y; c Y; modified 

cc <- 0.90
dd <- sqrt ((1+cc)/(1-cc)) # 4.359

thetaval <- seq(-8,8,by=0.01)
risk1val <- 1 + 0*thetaval # simple classic Y
risk2val <- cc^2 + (1-cc^2)*thetaval^2 # the estimator c Y

thetahat <- function(y)
{
ok <- 1*(abs(y) <= dd)
cc*y*(ok==1) + y*(ok==0)
}

curve(thetahat,-7,7)

risk3val <- 0*thetaval
for (j in 1:length(thetaval))
{
thetanow <- thetaval[j]
#
inte <- function(y)
{ (thetahat(y) - thetanow)^2 * dnorm(y-thetanow) }
#
aux1 <- integrate(inte,-100,0)$value
aux2 <- integrate(inte,0,100)$value
risk3val[j] <- aux1 + aux2 
}

par(bty="l")
par(cex.lab=1.4142) 

matplot(thetaval,cbind(risk1val,risk2val,risk3val),
   type="l",lwd=3,ylim=c(0.50,2),
   xlab="theta",ylab="three risk functions")