x = 1:length(wooster$data[,2])
usin = function(x, a, b, d)
{
        a + b * sin(((x - d) * 2 * pi)/365.25)
}
wu = usin(x, -30, 25, -75)
winter = c(rep(c(rep(1, 61), rep(0, 273), rep(1, 31)), 5),
        1)
spring = c(rep(c(rep(0, 61), rep(1, 91), rep(0, 365 - 91 -
        61)), 5), 0)
summer = c(rep(c(rep(0, 61 + 91), rep(1, 91), rep(0, 365 -
        91 - 61 - 91)), 5), 0)
fall = c(rep(c(rep(0, 61 + 91 + 91), rep(1, 91), rep(0, 365 -
        91 - 61 - 91 - 91)), 5), 0)
rescale.covariate  = function(x)
{
r.x = range(x)
x.01 = (x-r.x[1])/diff(r.x)
2*x.01-1
}
ydat = cbind(sin((x * 2 * pi)/365.25), cos((x * 2 * pi)/365.25
        ), rescale.covariate(x), winter, spring, summer, fall)

fit.model1 = pp.fit( - wooster$data[,2], threshold = wu,muinit=20.67317, siginit=13.05087 ,shinit=-0.0765431)

fit.model2 = pp.fit( - wooster$data[,2], threshold = wu, ydat = ydat, mul = 1:2,siglink = exp, muinit=c(-12,7,24), siginit=3,shinit= -0.1)

fit.model3 = pp.fit( - wooster$data[,2], threshold = wu, ydat = ydat,mul = 1:2, sigl = 1:2, siglink = exp,muinit=c(-15,9,29),siginit=c(0.4,0.3,0.1),shinit=-0.3)

fit.model4 = pp.fit( - wooster$data[,2], threshold = wu, ydat = ydat,mul = 1:2, shl=1:2,sigl = 1:2, siglink = exp,muinit=c(-15,9,29),siginit=c(0.4,0.3,0.1),shinit=c(0.1,0.1,0.1))