#### Simuliere Daten
n <- 1000
b <- 600
u <- rnorm(n)
u[(b+1):n] <- u[(b+1):n]*1.2
y <- arima.sim(model=list(ar=0.87,ma=-0.2),n=n,innov=u)

plot.ts(y)


#### Likelihood Funktion

LogLik <- function(theta, x, P, Q, b, restrict=TRUE){
n <- length(x)
n1 <- b
n2 <- n-b

c <- theta[1]
ar <- theta[2:(P+1)]
ma <- theta[(P+2):(P+Q+1)]
mu <- c/(1-sum(ar))

x <- x - mu

A <- matrix(0,P,P)
A[1,] <- ar
if (P>=2) A[2:P,1:(P-1)] <- diag(P-1)
if(max(abs(eigen(A)$values))>=1) return(Inf)

eta <-  na.omit(filter(x,c(1,-ar),sides=1))
u <- na.omit(filter(eta,c(-ma),"recursive"))
n <- length(u)

if (!restrict){
sigma1_sq <- theta[P+Q+2]
sigma2_sq <- theta[P+Q+3]
if (any(c(sigma1_sq,sigma2_sq)<=0)) return(Inf)
LL <- 0.5*(n1*log(sigma1_sq)+n2*log(sigma2_sq)+
  sum(u[1:n1]^2/sigma1_sq)+sum(u[(n1+1):n]^2/sigma2_sq))
} else {
sigma_sq <- theta[P+Q+2]
if (sigma_sq<=0) return(Inf)
LL <- 0.5*(n*log(sigma_sq)+sum(u^2/sigma_sq))
}

return(LL)
}






#est1 <-optim(f=LogLik, p=theta1, x=y, P=3, Q=3, b=600, control=list(maxit=5000),
#  restrict=FALSE)
#est1
#theta2 <- c(0,0.8,0.3,1)
#est2 <- optim(f=LogLik, p=theta2, x=y, P=1, Q=1, b=600, restrict=TRUE)
#est2
#LR <- 2*(est2$value-est1$value)
#pval <- pchisq(LR,df=1,lower.tail=FALSE)
#pval

P <- 3
Q <- 3
theta1 <- c(0,0.2,0.1,0.1,0.1,0.1,0.3,1,1)
est1 <- nlm(f=LogLik, p=theta1, x=y, P=3, Q=3, b=600,
  restrict=FALSE, iterlim = 1000, hessian=TRUE)
est1

estimates <- est1$estimate
stderrors <- sqrt(diag(solve(est1$hessian)))

result <- cbind(estimates,stderrors)
rownames(result) <- c("c",paste("ar",1:P),paste("ma",1:Q),"Var1","Var2")
result

theta2 <- c(0,0.8,0.3,1)
est2 <- nlm(f=LogLik, p=theta2, x=y, P=1, Q=1, b=600, restrict=TRUE)
est2

LR <- 2*(est2$minimum-est1$minimum)
pval <- pchisq(LR,df=1,lower.tail=FALSE)
pval



### Likelihood mit mehreren Bruchpunkten

LogLik <- function(theta, x, P, Q, b){


c <- theta[1]
ar <- theta[2:(P+1)]
ma <- theta[(P+2):(P+Q+1)]
mu <- c/(1-sum(ar))
x <- x - mu

A <- matrix(0,P,P)
A[1,] <- ar
if (P>=2) A[2:P,1:(P-1)] <- diag(P-1)
if(max(abs(eigen(A)$values))>=1) return(Inf)

eta <-  na.omit(filter(x,c(1,-ar),sides=1))
u <- na.omit(filter(eta,c(-ma),"recursive"))

b <- c(0,b,length(u))
s <- length(b)

sigma_sq <- theta[-(1:(P+Q+1))]
if (any(sigma_sq<=0)) return(Inf)

LL <- 0
for (i in 2:s) LL <- LL + 0.5*((b[i]-b[i-1])*log(sigma_sq[i-1])+
  sum(u[(b[i-1]+1):b[i]]^2/sigma_sq[i-1]))
return(LL)
}


b <- NULL
P <- 1
Q <- 1
par0 <- c(0.1,0.5,0.3,2)
e1 <- nlm(LogLik,p=par0,x=y,b=b,P=P,Q=Q)


LR <- 2*(est2$minimum-est1$minimum)
pval <- pchisq(LR,df=4,lower.tail=FALSE)
pval


qchisq(0.95,df=4)