R Under development (unstable) (2024-02-17 r85935 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> require("dse")
Loading required package: dse
Loading required package: tfplot
Loading required package: tframe

Attaching package: 'dse'

The following objects are masked from 'package:stats':

    acf, simulate

> Sys.info()
       sysname        release        version       nodename        machine 
     "Windows"   "Server x64"  "build 20348"     "CRANWIN3"       "x86-64" 
         login           user effective_user 
        "CRAN"         "CRAN"         "CRAN" 
> DSEversion()
       setRNG        tframe           dse 
   "2024.2-1" "2015.12-1.1"    "2024.2-1" 
>  
> fuzz <- 1e-6
> digits <- 18
> all.ok <- TRUE  
> 
> test.rng <- list(kind="Wichmann-Hill",seed=c(979,1479,1542),normal.kind="Box-Muller")
> 
> 
> # examples from Olivier Briet
> 
> 
> #fix constants:
> arma.model <- ARMA(A=c(1, -0.8), B=c(1,0.5), C=NULL, TREND=NULL)
> simulated.data <- simulate(arma.model, sampleT=1000, rng=test.rng)
> proposed.arma.model <- ARMA(A=c(1, -0.1), B=c(1,0.2), C=NULL, TREND=NULL,
+        constants=list(B=array(c(TRUE,TRUE),c(2,1,1))))
> 
> estMaxLik.result <-estMaxLik(simulated.data, proposed.arma.model)
> 
> print(estMaxLik.result)
neg. log likelihood= 1439.509 

A(L) =
1-0.8576495L1    

B(L) =
1+0.2L1    
> print(estMaxLik.result$estimates$like[1], digits=18)
[1] 1439.50913062493464
> print(TSmodel(estMaxLik.result)$B, digits=18)
, , 1

                     [,1]
[1,] 1.000000000000000000
[2,] 0.200000000000000011

> 
> good <-  c(1439.50913062493419, 1, 0.2) 
> error <- max(abs(good - 
+     c(estMaxLik.result$estimates$like[1],TSmodel(estMaxLik.result)$B)))
> cat("max. error ", max(error), "\n")
max. error  4.547474e-13 
>  
> if (any(is.na(error)) || any(is.nan(error)) || fuzz < error) all.ok <- FALSE  
>  
> 
> # trend (p-vector form):
> arma.model <- ARMA(A=c(1, -0.8), B=c(1,0.5), C=NULL, TREND=array(c(10),c(1,1)))
> simulated.data <- simulate(arma.model, sampleT=1000, rng=test.rng)
> estMaxLik.result <- estMaxLik(simulated.data, arma.model)
> 
> print(estMaxLik.result)
neg. log likelihood= 1403.508 
TREND=  10.70319
A(L) =
1-0.7864813L1    

B(L) =
1+0.4838096L1    
> print(estMaxLik.result$estimates$like[1], digits=18)
[1] 1403.50786987187303
> print(TSmodel(estMaxLik.result)$TREND, digits=18)
[1] 10.7031874021288669
> 
> good <-  c(1403.50786987187257 ,10.7031874021547484)
> error <- max(abs(good - 
+     c(estMaxLik.result$estimates$like[1], TSmodel(estMaxLik.result)$TREND)))
> cat("max. error ", max(error), "\n")
max. error  2.588152e-11 
>  
> if (any(is.na(error)) || any(is.nan(error)) || fuzz < error) all.ok <- FALSE  
>  
> 
>  
> # trend (matrix form):
> 
> arma.model<- ARMA(A=c(1, -0.8), B=c(1), C=NULL, TREND=matrix(c(1:10),10,1))
> 
> #N.B.: The estMaxLik computing time with matrix trend is very long for longer series
> simulated.data <- simulate(arma.model, sampleT=10, rng=test.rng)
> 
> estMaxLik.result <- estMaxLik(simulated.data, arma.model)
> 
> print(estMaxLik.result)
neg. log likelihood= 1.57484 
TREND=   1.000000  4.120202  3.738710  3.731154  4.361540  6.049436  8.028798  8.293034  6.691231 11.472080
A(L) =
1-0.818542L1    

B(L) =
1    
> print(estMaxLik.result$estimates$like[1], digits=18)
[1] 1.57484020410360159
> print(TSmodel(estMaxLik.result)$TREND, digits=18)
 [1]  1.00000000000000000  4.12020240437844976  3.73871024690458054
 [4]  3.73115357754047849  4.36154038371968511  6.04943599450627989
 [7]  8.02879756914381559  8.29303430382134898  6.69123127512212612
[10] 11.47207950374356500
> 
> good <- c( 1.57484020410352699,
+     1.0,  4.12020240399457460, 3.73871024870632374,  3.73115357779961387,
+           4.36154038369332131, 6.04943599227131035,  8.02879757311711373,
+           8.29303430238073780, 6.69123127611127888, 11.47207950229850937)
> 
> error <- max(abs(good -
+   c( estMaxLik.result$estimates$like[1],TSmodel(estMaxLik.result)$TREND)))
> 
> cat("max. error ", max(error), "\n")
max. error  3.973298e-09 
>  
> if (any(is.na(error)) || any(is.nan(error)) || fuzz < error) all.ok <- FALSE  
>  
>  
> if (! all.ok) stop("some tests FAILED")
> 
>  
> 
> proc.time()
   user  system elapsed 
   0.78    0.10    0.85