tests/TruePar.R

#  Tests of calculated "true values" for simulations.
#  See notes in ~/papers/InProgress/TimeSeriesFA/calculations.R

require("tsfa")

if (require("CDNmoney")) {

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

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

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)




##poppar

DX <- diff(MBcomponents)
stds <- sqrt(diag(cov(DX)))

poppar <- factanal(factors = 2, covmat = cov(DX), n.obs = Tobs(DX),
            scores = "none", rotation = "none")

# Quartimin (= oblimin with parameter 0) rotation with Kaiser normalization
rownorms  <- sqrt(rowSums(poppar$loadings^2))

KBorth    <- diag(1/rownorms) %*% poppar$loadings
rotpoppar <- GPFoblq(A = KBorth, Tmat = diag(1, 2), method="quartimin")
Boblq     <- diag(stds * rownorms) %*% rotpoppar$loadings
PhiOblq   <- rotpoppar$Phi
Omega     <- diag(stds * poppar$uniquenesses * stds)
Psi       <- 0.5 * Omega

# So the true parameters in the simulation are
  Boblq
  PhiOblq
  Omega
  Psi

SigmaOblq <- Boblq %*% PhiOblq %*% t(Boblq) + Omega

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))

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


###  Tests to check that calculated values have not changed

#  NB  fuzz is relaxed even more further below
fuzz <- 1e-5 # this is pretty relaxed compared to my usual tests, eventually
#              tighten up with values printed to higher precision.

# ######################################################################
# # How well do the means fit?
 
 meanDX  <- colMeans(diff(MBcomponents))
 SigmaOblq <- Boblq %*% PhiOblq %*% t(Boblq) + Omega

# kappaOblq <- solve(t(Boblq) %*% solve(SigmaOblq) %*% Boblq) %*%
#		      t(Boblq) %*% solve(SigmaOblq) %*% t(meanDX)
 kappaOblq <- solve(t(Boblq) %*% solve(SigmaOblq, Boblq),
		      t(Boblq)) %*% solve(SigmaOblq) %*% meanDX
 
 meanhat <- t(Boblq %*% kappaOblq)
 

######################################################################

 tst <- c(  2.882879, -0.3320896, 0.1614583, 10.56027, 9.028839, 14.88725)
 if( fuzz < max(abs(meanhat - tst))) {
    cat("Calculated value is not the same as test value in test 1. Value:\n")
    printTestValue(meanhat, digits=18)
    cat("difference:\n")
    print(meanhat - tst, digits=18)
    all.ok <- FALSE  
    }
    
 tst <- c(  1.623936, 0.8437919, 2.553238, 10.21236, 9.243064, 56.72988)
 if( fuzz < max(abs( meanDX - tst ))) {
    cat("Calculated value is not the same as test value in test 2. Value:\n")
    printTestValue(meanDX, digits=18)
    cat("difference:\n")
    print(meanDX - tst, digits=18)
    all.ok <- FALSE  
    }


 tst <- c(  77.524138564771, -139.35681520788, -93.676330196783, 3.4067233632789,
      -2.3176922549429, -73.757662639795)
 if( fuzz < max(abs( 100*(meanhat-meanDX)/meanDX - tst ))) {
    cat("Calculated value is not the same as test value in test 3. Value:\n")
    printTestValue(100*(meanhat-meanDX)/meanDX, digits=18)
    cat("difference:\n")
    print(100*(meanhat-meanDX)/meanDX - tst, digits=18)
    all.ok <- FALSE  
    }

 
              
# Check correctness
DetaTrue <- diff(etaTrue, lag=1)
SigDeta  <- cov(DetaTrue)
CheckFax <- 100 * (SigDeta - PhiOblq)/PhiOblq   # relative difference in %

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



 tst <- c(  0.40629016392777, 0.63131919291667, 0.74051211937755, 0.27589064829569,
         0.80375582540918, 0.84226727885335)
 if( fuzz < max(abs(  poppar$uniquenesses - tst))) {
    cat("Calculated value is not the same as test value in test 5. Value:\n")
    printTestValue(poppar$uniquenesses, digits=18)
    cat("difference:\n")
    print(poppar$uniquenesses - tst, digits=18)
    all.ok <- FALSE  
    }


 tst <- t(matrix(c( 
       0.76269855837340,  0.109545890987616,
       0.38940769411444, -0.465876126560372,
       0.37235218322458, -0.347621664499332,
       0.85015261178239,  0.036737207574017,
       0.13286428348581,  0.422603915908731,
       0.17306303686227,  0.357462518428758), 2,6))
 if( fuzz < max(abs(  poppar$loadings - tst ))) {
    cat("Calculated value is not the same as test value in test 6. Value:\n")
    printTestValue(poppar$loadings[,], digits=18)
    cat("difference:\n")
    print(poppar$loadings[,] - tst, digits=18)
    all.ok <- FALSE  
    }


#poppar
#		Factor1 Factor2
#SS loadings	  1.642   0.658
#Proportion Var   0.274   0.110
#Cumulative Var   0.274   0.383

fuzz <- 1e-6 

## simulation values

#tst <- diag(c(72.63349, 1233.026, 87.33772, 629.3927, 3967.983, 12163.26))
tst <- diag(c(72.633490218431234, 1233.026245431895177, 87.337721037020572,
              629.392699084312198, 3967.982989812266169, 12163.258995566555313))
 if( fuzz < max(abs( Omega - tst ))) {
    cat("Calculated value is not the same as test value in test 7. Value:\n")
    printTestValue(diag(Omega), digits=18)
    cat("difference:\n")
    print(Omega - tst, digits=18)
    all.ok <- FALSE  
    }


# using early version of GPA 
# tst <- t(matrix(c( 
#	8.8424730937083,   5.2034756517999,
#      23.8239001815301, -12.5767329305997,
#	5.1878379563914,  -1.9710232106657,
#      36.7834441658135,  16.9430523091825,
#      -2.8494840729111,  31.0224249106556,
#	2.6047424392098,  47.6304266899614), 2,6))

 tst <- t(matrix(c( 
      8.84252375520804534,   5.20339025247098252, 
     23.8237818511178361,  -12.576960598600591, 
      5.18781954415890123,  -1.97107285697162915, 
     36.7836098529645881,   16.9426974878479477, 
     -2.8491870985089478,   31.0224494555067913, 
      2.60519903617723481,  47.6303973396630411), 2,6))

 if( fuzz < max(abs( Boblq - tst ))) {
    cat("Calculated value is not the same as test value in test 7.5. Value:\n")
    printTestValue(Boblq, digits=18)
    cat("difference:\n")
    print(Boblq - 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")

Try the tsfa package in your browser

Any scripts or data that you put into this service are public.

tsfa documentation built on May 29, 2017, 6 p.m.