蒙特卡洛数据生成
N=1500#迭代次数
> burnin=500      #定义预烧期
> k<-0.15
> kappa<-0.003
> rho0<- -0.1
> rho1<- -0.3
> sigma0<-1.0
> sigma1<- 0.5
> beta0<- -0.5
> beta1<- -1.0
> alpha0<- 0.6
> alpha1<- 0.9
> x0<--0.0065
> x<-vector(mode="numeric",length=N)
> h<-vector(mode="numeric",length=N+1)
> epsilon<-c(rnorm(N,0,1))     # 定义为一行向量，每个元素都服从均值为0方差为1的正态分布
> lambda<-c(rbinom(N,1,kappa))   # 定义为一行向量，每个元素都是伯努利随机变量
> eta<-vector(mode="numeric",length=N+1)
> if (x0<=0)   h[1]<-rnorm(1,(beta0)/(1-(alpha0)),sqrt(1/(1-(alpha0)^2))) else
+ h[1]<-rnorm(1,(beta1)/(1-(alpha1)),sqrt(1/(1-(alpha1)^2)))
> x[1]<-k*lambda[1]+exp(h[1]/2)*epsilon[1]
> for (i in 2:N)
+ {
+ if (x[i-1]<=0) {
+ eta<-rnorm(1,(rho0)*(sigma0)*epsilon[i-1],(sigma0)*sqrt(1-(rho0)^2))
+ h<-beta0+alpha0*h[i-1]+eta
+ }  else {
+ eta<-rnorm(1,(rho1)*(sigma1)*epsilon[i-1],(sigma1)*sqrt(1-(rho1)^2))
+ h<-beta1+alpha1*h[i-1]+eta}
+ x<-k*lambda+exp(h/2)*epsilon
+ }
> if (x[N]<=0) {
+ eta[N+1]<-rnorm(1,(rho0)*(sigma0)*epsilon[N],(sigma0)*sqrt(1-(rho0)^2))
+ h[N+1]<-beta0+alpha0*h[N]+eta[N+1]
+ }  else {
+ eta[N+1]<-rnorm(1,(rho1)*(sigma1)*epsilon[N],(sigma1)*sqrt(1-(rho1)^2))
+ h[N+1]<-beta1+alpha1*h[N]+eta[N+1]}
> b<-burnin+1
> X<-x[b:N]
> H<-h[b:(N+1)]
构建似然函数
> LogL<-function(theta)
+ {beta0<-theta[1]
+ beta1<-theta[2]
+ alpha0<-theta[3]
+ alpha1<-theta[4]
+ sigma0<-theta[5]
+ sigma1<-theta[6]
+ rho0<-theta[7]
+ rho1<-theta[8]
+ k<-theta[9]
+ kappa<-theta[10]
+ n=length(X)
+ X0<-x[burnin]
+ epsilon<-c(rnorm(n,0,1))
+ if (X0<=0){
+ L1<-(kappa*(1/sqrt(2*pi*exp(H[1])))*exp((-(X[1]-k)^2)/(2*exp(H[1])))+(1-kappa)*(1/sqrt(2*pi*exp(H[1])))*exp((-(X[1])^2)/(2*exp(H[1]))))*(1/(sqrt(2*pi*(1-(rho0)^2))*sigma0))*exp((-(H[2]-beta0-alpha0*H[1]-sigma0*rho0*epsilon[1])^2)/(2*((sigma0)^2)*(1-(rho0)^2)))
+ } else {
+ L1<-(kappa*(1/sqrt(2*pi*exp(H[1])))*exp((-(X[1]-k)^2)/(2*exp(H[1])))+(1-kappa)*(1/sqrt(2*pi*exp(H[1])))*exp((-(X[1])^2)/(2*exp(H[1]))))*(1/(sqrt(2*pi*(1-(rho1)^2))*sigma0))*exp((-(H[2]-beta1-alpha1*H[1]-sigma1*rho1*epsilon[1])^2)/(2*((sigma1)^2)*(1-(rho1)^2)))
+ }
+ for(i in 2:n){
+ if (X[i-1]<=0){
+ Li<-kappa*(1/sqrt(2*pi*exp(H)))*exp((-(X-k)^2)/(2*exp(H)))+(1-kappa)*(1/sqrt(2*pi*exp(H)))*exp((-(X)^2)/(2*exp(H)))*(1/(sqrt(2*pi*(1-(rho0)^2))*sigma0))*exp((-(H[i+1]-beta0-alpha0*H-sigma0*rho0*epsilon)^2)/(2*((sigma0)^2)*(1-(rho0)^2)))
+ } else {
+ Li<-kappa*(1/sqrt(2*pi*exp(H)))*exp((-(X-k)^2)/(2*exp(H)))+(1-kappa)*(1/sqrt(2*pi*exp(H)))*exp((-(X)^2)/(2*exp(H)))*(1/(sqrt(2*pi*(1-(rho1)^2))*sigma1))*exp((-(H[i+1]-beta1-alpha1*H-sigma1*rho1*epsilon)^2)/(2*((sigma1)^2)*(1-(rho1)^2)))
+ }
+ }
+ ll<-sum(log(Li))
+ return( ll)
+ }
> result<-maxLik(LogL,start=c(-0.5,-1.0,-0.6,-0.9,1,0.5,-0.1,-0.3,3,0.005))
> print(result)
Maximum Likelihood estimation
Newton-Raphson maximisation, 1 iterations
Return code 3: Last step could not find a value above the current.
Boundary of parameter space?  这里是什么意思？
Consider switching to a more robust optimisation method temporarily.
Log-Likelihood: -2.223315 (10 free parameter(s))
Estimate(s): -0.5000001 -1 -0.6 -0.8999999 1 0.4999995 -0.09999986 -0.2999997 3 0.004999788
> 我估计出来的似然值也不对啊，找不到错误，求大神指教。