if (FALSE) {
## simulate the times of a Poisson process on [0,1] with intensity
## function nu(t)=100*(2+cos(2*pi*t)).
tms <- simPois(int=function(x)100*(2+cos(2*pi*x)),int.M=300)
## calculate a nonparametric estimate of the intensity function
int <- lpint::lpint(tms,Tau=1)
matplot(int$x,int$y+qnorm(0.975)*outer(int$se,c(-1,0,1)),type="l",lty=c(2,1,2),
col=1,xlab="t",ylab="nu(t)")
curve(100*(2+cos(2*pi*x)),add=TRUE,col="red")
## simulate an IHSEP on [0,1] with baseline intensity function
## nu(t)=100*(2+cos(2*pi*t)) and excitation function
## g(t)=0.5*8*exp(-8*t)
asep <- simHawkes1(nu=function(x)200*(2+cos(2*pi*x)),nuM=600,
g=function(x)8*exp(-16*x),gM=8)
## get the birth times of all generations and sort in ascending order
tms <- sort(unlist(asep))
## calculate the minus loglikelihood of an SEPP with the true parameters
mloglik1a(tms,TT=1,nu=function(x)200*(2+cos(2*pi*x)),
g=function(x)8*exp(-16*x),Ig=function(x)8/16*(1-exp(-18*x)))
## calculate the MLE for the parameter assuming known parametric forms
## of the baseline intensity and excitation functions
est <- optim(c(400,200,2*pi,8,16),
function(p){
mloglik1a(jtms=tms,TT=1,
nu=function(x)p[1]+p[2]*cos(p[3]*x),
g=function(x)p[4]*exp(-p[5]*x),
Ig=function(x)p[4]/p[5]*(1-exp(-p[5]*x)))
},
hessian=TRUE,control=list(maxit=5000,trace=TRUE))
## point estimate by MLE
est$par
## standard error estimates:
diag(solve(est$hessian))^0.5
}
Run the code above in your browser using DataLab