# Cut-and-Paste Code below into Window Above and RUN
require(dse)
require(matlab)
#
# Measurement Matrix
# CR1=(Overall Growth), CR2=(LU+KOF+EF -EG-Q-CO2), CR3=(EF-LU)
#
# EN.ATM.CO2E.KT NY.GDP.MKTP.KD EG.USE.COMM.KT.OE SL.TLF.TOTL.IN
#[1,] 0.3463 0.342 0.331 0.3524
#[2,] -0.2176 -0.296 -0.371 -0.1381
#[3,] 0.0605 -0.172 -0.299 -0.0653
# SP.POP.TOTL SL.UEM.TOTL.ZS KOF EF HDI
#[1,] 0.3523 0.252 0.339 0.320 0.3537
#[2,] -0.1083 0.715 0.231 0.363 0.0314
#[3,] -0.0113 -0.634 0.296 0.615 0.0757
#
# Fraction of Variance
#[1] 0.876 0.967 0.989 0.994 0.997 0.999 0.999 1.000 1.000
#
AIC <- function(model) {informationTestsCalculations(model)[3]}
f <- matrix( c( 0.98792655, 0.02900433, 0.1716011, 0.168154718,
-0.04332413, 1.02393353, 0.2438378, -0.031433144,
-0.01630805, 0.02497518, 0.8391857, -0.002772357,
0.000000000, 0.000000000, 0.000000, 1.00000
),byrow=TRUE,nrow=4,ncol=4)
#
# To stabilize the system, uncomment next line
# f[2,2] <- .9
#
# To turn Growth into a random walk with drift
# f[1,1] <- 1
#
h <- eye(3,4)
k <- f[1:4,1:3,drop=FALSE]
CRL20 <- SS(F=f,H=h,K=k,z0=c( 0.168154718, -0.031433144, -0.002772357, 1.0000000),
output.names=c("CR1","CR2","CR3"))
stability(SS(F=f[1:3,1:3,drop=FALSE],H=eye(3),Q=eye(3),R=eye(3)))
shockDecomposition(toSSChol(CRL20))
#tfplot(simulate(CRL20,sampleT=50,noise=matrix(0,50,2),start=1960))
CRL20.data <- simulate(CRL20,sampleT=50,start=1960)
m <- l(CRL20,CRL20.data)
#tfplot(m)
CRL20.f <- forecast(m,horizon=50)
tfplot(CRL20.f)
AIC(m)