tests/dse1tst07.R

In dse: Dynamic Systems Estimation (Time Series Package)

 require("dse")
 Sys.info()
 DSEversion()
 
fuzz.small <- 1e-14
fuzz.large <- 1e-10
digits <- 18
all.ok <- TRUE  


cat("dse test 7a...\n")
  d <-20
  true.roots <- c(-1/seq(d),1/seq(d),-seq(d),seq(d)) 
  A <- array(0, c(2,length(true.roots),length(true.roots)))
  A[1,,] <- diag(1,length(true.roots))
  A[2,,] <- diag(-true.roots, length(true.roots))
  # the following relies on roots using by.poly=FALSE
  # Splus needed options(expressions=1024)
   tst <-  sort(roots( ARMA(A=A, B=diag(1,length(true.roots)) ),by.poly=FALSE))
   good <- sort(true.roots)
   error <- max(Mod(good - tst))
   cat("max. error ", max(error))

   if (any(is.na(error)) || any(is.nan(error)) || fuzz.small < error ) 
     {printTestValue(c(tst), digits=18)
      all.ok <- FALSE  
     }

  #  complex roots test
cat("dse test 7b...\n")

  i <- pi*(1:10)/10.1  # this is half circle, but conjs also get generated
  # div by 10.1 instead of 10 prevents a real root with = conj
  true.roots <- complex(real=cos(i), imaginary=sin(i))  # on unit circle
  # scale simplifies sorting
  true.roots <- c(true.roots*(1+.2*1:10), true.roots*(1:10)/10) 
  A <- array(0, c(3,length(true.roots),length(true.roots)))
  A[1,,] <- diag(1,length(true.roots))
  A[2,,] <- diag(-2*Re(true.roots), length(true.roots))
  A[3,,] <- diag(Re(true.roots*Conj(true.roots)), length(true.roots))
  est.roots <- roots( ARMA(A=A, B=diag(1,length(true.roots)) ))
  ec <- 0<=Im(est.roots)

  z  <- est.roots[ ec][order(Mod(est.roots[ ec]))]
  zz <-  true.roots[order(Mod(true.roots))]

  error <- max(Mod(z - zz))
   cat("max. error ", max(error))

   if (any(is.na(error)) || any(is.nan(error)) || fuzz.small < error ) 
     {printTestValue(c(z), digits=18)
      printTestValue(c(zz), digits=18)
      all.ok <- FALSE  
     }

cat("dse test 7c...\n")

  z  <- est.roots[!ec][order(Mod(est.roots[!ec]))]
  zz <- Conj(true.roots)[order(Mod(true.roots))]
  error <- max(Mod(z - zz))
   cat("max. error ", max(error))

   if (any(is.na(error)) || any(is.nan(error)) || fuzz.small < error ) 
     {printTestValue(c(z), digits=18)
      printTestValue(c(zz), digits=18)
      all.ok <- FALSE  
     }


  if (! all.ok) stop("some tests FAILED")