tests/MCestTSF_ML.R

In tsfa: Time Series Factor Analysis

#  See notes in file calculations.R for TimeSeriesFA paper 

require("tsfa")

if (require("CDNmoney")) {

data("CanadianMoneyData.asof.6Feb2004", package="CDNmoney")

# for monte carlo
require("dse")
require("EvalEst")  # for EstEval

fuzz <- 1e-8
all.ok <- TRUE  

### Construct data

cpi <- 100 * M1total / M1real
seriesNames(cpi) <- "CPI"
popm <- M1total / M1PerCapita
seriesNames(popm) <- "Population of Canada"

z <- tframed(tbind(
    MB2001,
    MB486 + MB452 + MB453 ,
    NonbankCheq,
    MB472 + MB473 + MB487p,
    MB475,
    NonbankNonCheq + MB454 + NonbankTerm + MB2046 + MB2047 + MB2048 +
    MB2057 + MB2058 + MB482),
    names=c("currency", "personal cheq.", "NonbankCheq",
    "N-P demand & notice", "N-P term", "Investment")
    )

z <- tfwindow(z, start=c(1986,1))
if( all(c(2003,12) ==end(z))) z <-tfwindow(z, end=c(2003,11))
MBcomponents <- 1e8 * z/matrix(tfwindow(popm * cpi,tf=tframe(z)),Tobs(z),6)

### Specify "true" parameters and factors

Omega <- diag(c(72.633490218431234, 1233.026245431895177, 87.337721037020572,
              629.392699084312198, 3967.982989812266169, 12163.258995566555313))

Boblq <- t(matrix(c( 
     8.8424730199723260,   5.2034757439511159,
    23.8239003553122046, -12.5767326858555819,
     5.1878379834837549,  -1.9710231572687940,
    36.7834439249370746,  16.9430526918934632,
    -2.8494845070603847,  31.0224248853059343,
     2.6047417719514878,  47.6304267232332990), 2,6))

PhiOblq <- t(matrix(c(
    1.0000000000000002220, 0.0094910545788177599,
    0.0094910545788177599, 1.0000000000000002220),2,2))

etaBart <- MBcomponents %*% solve(Omega) %*% Boblq %*% (
              solve( t(Boblq) %*% solve(Omega) %*% Boblq ) )

DetaBart <- diff(etaBart, lag=1)
SDE      <- cov(DetaBart)       
RR1 <- chol(SDE)      # upper triangular: SDE = RR1' RR1
RR2 <- chol(PhiOblq)  # ditto
PP  <- t(RR2) %*% solve(t(RR1))
Psi       <- 0.5 * Omega

etaTrue <- tframed(etaBart %*% t(PP), tf=tframe(MBcomponents))


rngValue10 <- list(seed=10, kind="Mersenne-Twister", normal.kind="Inversion")

cat("Starting the trucated monte carlo experiment test...\n")

EE.ML5 <- EstEval(TSFmodel(Boblq, f=etaTrue, positive.measures=FALSE),
    replications=5, rng=rngValue10, quiet=FALSE,
    simulation.args=list(Cov=Psi),
    estimation="estTSF.ML", 
    estimation.args=list(2, BpermuteTarget=Boblq),
    criterion ="TSFmodel")

mc5SummaryML <- summary(EE.ML5)

# the next were not tested with early version of GPA  

 tst <- c(-6.35811456287209165,  0.0607392955938533957)
#  above was with noIC=TRUE
 tst <- c(-2.92588953587939571,  -4.71404513263362457)

 if( fuzz < max(abs( (mc5SummaryML$meanhatf.error)  - tst ))) {
    cat("Calculated value is not the same as test value in test 19. Value:\n")
    printTestValue(mc5SummaryML$meanhatf.error, digits=18)
    cat("difference:\n")
    print(mc5SummaryML$meanhatf.error - tst, digits=18)
    all.ok <- FALSE  
    }

 tst <-  c(11.1176063210974636,  18.1248814664440161)
#  above was with noIC=TRUE
 tst <-  c(8.60111941828624538,  24.7872684764427618)

 if( fuzz < max(abs( mc5SummaryML$meanSDhatf  - tst ))) {
    cat("Calculated value is not the same as test value in test 21. Value:\n")
    printTestValue(mc5SummaryML$meanSDhatf, digits=18)
    cat("difference:\n")
    print(mc5SummaryML$meanSDhatf - tst, digits=18)
    all.ok <- FALSE  
    }

 tst <- t(matrix(c( 
      -1.07625157111048164 , 0.703342731082228667, 
      -1.2178110662245345  , 4.70567370637503757 , 
      -0.167039105282395184, 0.226947650291669634, 
      -6.57680144713607007 , 6.63837150437990431 , 
      -7.64998746817823871 , -1.11899415928688128 , 
      -4.34709518945724671 , 3.75421190419606887 ), 2,6))
#  above was with noIC=TRUE
 tst <- t(matrix(c( 
      -1.06570250190681648 , -0.0130131676780012029 , 
      -1.85077437881246709 , 3.15697064361194002 , 
      -0.421861915762349327 , -0.337332002481196769 , 
      -2.95530986474609847 , -0.00875023040947908726 , 
      1.05168217191284175 , 2.10624959737685558 , 
      2.07817788622315858 , -7.99355084897341328), 2,6))

 if( fuzz < max(abs( mc5SummaryML$meanhatB.error - tst ))) {
    cat("Calculated value is not the same as test value in test 23. Value:\n")
    printTestValue(mc5SummaryML$meanhatB.error, digits=18)
    cat("difference:\n")
    print(mc5SummaryML$meanhatB.error - tst, digits=18)
    all.ok <- FALSE  
    }

 tst <- t(matrix(c( 
      2.09171817126895432 , 3.10870214359533614 , 
      3.31424639812060606 , 9.86511334745889101 , 
      0.510677029668560656 , 1.93445199388070876 , 
      8.13415536978554243 , 11.7227922878270903 , 
     13.0294442536303503 , 6.96160177068301955 , 
     10.7630330480164904 , 19.6839038208302632), 2,6))
#  above was with noIC=TRUE
 tst <- t(matrix(c( 
      2.00786152578053922 , 3.51641487208566783 , 
      5.55688084114462377 , 9.53701794115509216 , 
       1.2949230299373764 , 1.46529585264552353 , 
      5.85468954626221016 , 10.3962361962235175 , 
        14.68586418448686 , 7.04445496431349572 , 
      11.7090617911295691 , 19.4653799970791646), 2,6))

 if( fuzz < max(abs( mc5SummaryML$SDhatB - tst ))) {
    cat("Calculated value is not the same as test value in test 24. Value:\n")
    printTestValue(mc5SummaryML$SDhatB, digits=18)
    cat("difference:\n")
    print(mc5SummaryML$SDhatB - tst, digits=18)
    all.ok <- FALSE  
    }


 fuzz <- 1e-5
 
 tst <- c(-0.0268800596792208475,  -0.0446045635023343962)
#  above was with noIC=TRUE
 tst <- c( -0.00762399780916277366,  -0.0437505013835693748)

 if( fuzz < max(abs(  (mc5SummaryML$meanhatDf.error)  - tst ))) {
    cat("Calculated value is not the same as test value in test 27. Value:\n")
    printTestValue(mc5SummaryML$meanhatDf.error, digits=18)
    cat("difference:\n")
    print(mc5SummaryML$meanhatDf.error - tst, digits=18)
    all.ok <- FALSE  
    }

 tst <- c( 0.689731516283344948,  0.80395876793842791)
#  above was with noIC=TRUE
 tst <- c(0.622493027392810316,  0.878995874209845596)

 if( fuzz < max(abs( mc5SummaryML$meanSDhatDf  - tst ))) {
    cat("Calculated value is not the same as test value in test 29. Value:\n")
    printTestValue(mc5SummaryML$meanSDhatDf, digits=18)
    cat("difference:\n")
    print(mc5SummaryML$meanSDhatDf - tst, digits=18)
    all.ok <- FALSE  
    }



#### this is fairly long 

cat("Starting the main monte carlo experiment...\n")

EE.ML100 <- EstEval(TSFmodel(Boblq, f=etaTrue, positive.measures=FALSE),
    replications=100, rng=rngValue10, quiet=FALSE,
    simulation.args=list(Cov=Psi),
    estimation="estTSF.ML", 
    estimation.args=list(2, BpermuteTarget=Boblq),
    criterion ="TSFmodel")

mc100SummaryML <- summary(EE.ML100)

fuzz <- 1e-5  # there seems to be a fair amount of "wobble" in these results 


# tst <- c(-1.437001,        -7.074382)
# tst <- c(-1.4366371949955, -7.0749007937085)
# using early version of GPA 
# tst <- c(-1.4366366998990,     -7.0749017751486)
  tst <- c(-1.43661535460228751, -7.07498653794059429)
#  above was with noIC=TRUE
  tst <- c(-4.25581051242354835,  -6.55708042849400741) # Linux Mandrake 32bit
  tst <- c(-4.25580918879642311,  -6.55708238950199629) # Solaris
  tst <- c(-4.2558098,            -6.5570813) 

 if( fuzz < max(abs( (mc100SummaryML$meanhatf.error)  - tst ))) {
    cat("Calculated value is not the same as test value in test 19. Value:\n")
    printTestValue(mc100SummaryML$meanhatf.error, digits=18)
    cat("difference:\n")
    print(mc100SummaryML$meanhatf.error - tst, digits=18)
    all.ok <- FALSE  
    }

# tst <- c(10.72659,        20.42743)
# using early version of GPA 
# tst <-  c(10.726408029255,     20.427522660190)
  tst <-  c(10.7266141051993671, 20.4274178936896966)
#  above was with noIC=TRUE
  tst <-  c(12.2212706117814953,  21.2410851409922579)

 if( fuzz < max(abs( mc100SummaryML$meanSDhatf  - tst ))) {
    cat("Calculated value is not the same as test value in test 21. Value:\n")
    printTestValue(mc100SummaryML$meanSDhatf, digits=18)
    cat("difference:\n")
    print(mc100SummaryML$meanSDhatf - tst, digits=18)
    all.ok <- FALSE  
    }

# using early version of GPA 
# tst <- t(matrix(c( 
#      -0.43072700321218, -0.669882007666349,
#      -2.45915806739033,  0.443541199032001,
#      -0.39598579113561, -0.055347774362779,
#      -1.28804986891355, -2.462726422671402,
#       1.03579550708264,  0.413078009465796,
#       2.50999408754555, -0.359689824353282), 2,6))

 tst <- t(matrix(c(
      -0.430727502635209092, -0.669890058809444078 , 
      -2.45916993900246084 ,  0.443526818436103198 , 
      -0.395990134050530962, -0.0553539782496834665, 
      -1.28804570353536718 , -2.46275799746584845  , 
       1.0358428846874721  , 0.413067269627077849  , 
       2.51004726054249838 , -0.359679178084249429), 2,6))
#  above was with noIC=TRUE
 tst <- t(matrix(c(
      -0.632879784699293069 , -0.354831582672403556 , 
      -1.59531543933821141 , 0.582831484242154474 , 
      -0.399060012837218814 , -0.0892344631765265017 , 
      -1.83055582409809148 , -0.331659929649745777 , 
      1.15013393059134073 , -1.14983030160259503 , 
      0.402743877294618802 , 1.06826205318063217 ), 2,6)) # Linux Mandrake 32bit
 tst <- t(matrix(c( 
      -0.632879702135179301 , -0.354831965108128422 , 
      -1.59531627367130113 , 0.582832449049474732 , 
      -0.399060005892669878 , -0.0892348968104166307 , 
      -1.83055540872476286 , -0.331661458302995982 , 
      1.15013525749168988 , -1.14983284599435365 , 
      0.402743928085540492 , 1.06826581647906949), 2,6))# Solaris
fuzz <- 1e-5  

 if( fuzz < max(abs( mc100SummaryML$meanhatB.error - tst ))) {
    cat("Calculated value is not the same as test value in test 23. Value:\n")
    printTestValue(mc100SummaryML$meanhatB.error, digits=18)
    cat("difference:\n")
    print(mc100SummaryML$meanhatB.error - tst, digits=18)
    all.ok <- FALSE  
    }

# using early version of GPA 
# tst <- t(matrix(c( 
#	 2.1177074120071,  2.8630417742169,
#	 5.8729782096725,  8.9365319798381,
#	 1.4920603854630,  1.9177841472202,
#	 6.9114403738557, 11.0123879149370,
#	 9.6088971180058, 10.4635264914877,
#	12.4426965380737, 17.4104141305891), 2,6))

 tst <- t(matrix(c( 
      2.11774533901723538,  2.86303620239022649, 
      5.87290229669339681,  8.9365652774522939 , 
      1.49204785812003604,  1.91779395295759936, 
      6.91158798109280426, 11.0123561946326376 , 
      9.60890954609334536, 10.463541675796586  , 
     12.4426951488466564,  17.4103921927887733 ),2,6))
#  above was with noIC=TRUE
 tst <- t(matrix(c( 
      2.33840665475836795 , 3.13852181903253502 , 
      6.66951161084711774 , 8.71483640509499935 , 
      1.33929845112549506 ,  2.0031133449915548 , 
      8.52324682951719304 , 11.2751832187060881 , 
      10.5144095556300563 , 7.63423805298109581 , 
      12.7681497781142195 , 15.1413132072165837 ),2,6))# Linux Mandrake 32bit
 tst <- t(matrix(c( 
      2.33840667415145909 , 3.13852213524006274 , 
      6.66951154478138974 ,  8.7148370261171344 , 
      1.33929845221169086 , 2.00311359935738809 , 
      8.52324720810181091 , 11.2751832991476277 , 
      10.5144106848696381 , 7.63424186177272812 , 
      12.7681497223306177 ,   15.14131310687468),2,6))# Solaris


 if( fuzz < max(abs( mc100SummaryML$SDhatB - tst ))) {
    cat("Calculated value is not the same as test value in test 24. Value:\n")
    printTestValue(mc100SummaryML$SDhatB, digits=18)
    cat("difference:\n")
    print(mc100SummaryML$SDhatB - tst, digits=18)
    all.ok <- FALSE  
    }


 fuzz <- 1e-5

# tst <-  c(0.01709371,        1.211911)
# tst <-  c(0.017095808203623, 1.1266841304661)
# above was when Df was percent change (now PCf) rather than difference

 tst <- c(0.0079738515556172, -0.043994019995693)
#  above was with noIC=TRUE
 tst <- c( -0.00505274415462879378,  -0.042964869565685973)

 if( fuzz < max(abs(  (mc100SummaryML$meanhatDf.error)  - tst ))) {
    cat("Calculated value is not the same as test value in test 27. Value:\n")
    printTestValue(mc100SummaryML$meanhatDf.error, digits=18)
    cat("difference:\n")
    print(mc100SummaryML$meanhatDf.error - tst, digits=18)
    all.ok <- FALSE  
    }

#  tst <- c(1.367774,        94.06914)
#  tst <- c(1.3677449004785, 92.129604373376)
# above was when Df was percent change (now PCf) rather than difference
 tst <- c( 0.67484628813627, 1.0077586888048)
#  above was with noIC=TRUE
 tst <- c(0.696639940140005298,  0.994209984498502064)

 if( fuzz < max(abs( mc100SummaryML$meanSDhatDf  - tst ))) {
    cat("Calculated value is not the same as test value in test 29. Value:\n")
    printTestValue(mc100SummaryML$meanSDhatDf, digits=18)
    cat("difference:\n")
    print(mc100SummaryML$meanSDhatDf - tst, digits=18)
    all.ok <- FALSE  
    }

cat("tests completed.\n")

if (! all.ok) stop("some tests FAILED")
} else cat("CDNmoney data not available. Tests skipped.\n")