NA_RE

#

# Cut-and-Paste code below into Window Above and Run

# NA_RE North America Roman Empire Model

#

# WRE World Model Roman Empire  Input (RE)

#

# Measurement Matrix 

#         Q      N      T

#[1,] 0.577  0.578 -0.577

#[2,] 0.528  0.274  0.804

#[3,] 0.622 -0.769 -0.147

#

# Fraction of Variance 

#[1] 0.998 1.000 1.000

#

merge.forecast <- function (fx,n=1) {

x <- splice(fx$pred,fx$forecast[[n]])

colnames(x) <- seriesNames(fx$data$output)

return(x)

}

AIC <- function(model) {informationTestsCalculations(model)[3]}

require(dse)

require(matlab)

f <- matrix( c(9.990858e-01, 6.700905e-03,  3.284895e+09,  1.173565e-02,

                    -6.563034e-04, 1.004810e+00,  2.358134e+09, -1.717386e-04,

                     8.324437e-19, 2.223337e-17, -7.995155e-02 , 1.087672e-15,

0.00000000,  0.0000000,  0.0000000,  1.0000000000

),byrow=TRUE,nrow=4,ncol=4)

h <- eye(3,4)

k <- (f[,1:3,drop=FALSE])

WRE <- SS(F=f,H=h,K=k,z0=c(1.173565e-02, -1.717386e-04 , 1.087672e-15,  1.0000000000),

              output.names=c("W1","W2","W3"))

WRE.data <- simulate(WRE,sampleT=150,noise=matrix(0,150,3))

WRE.f <- forecast(m <- l(WRE,WRE.data),horizon=150)

WRE.fx <- merge.forecast(WRE.f)

#

#

# NA_RE Latin Model Roman Empire  (RE)

#

# Measurement Matrix 

#         Q      N

#[1,] 0.707  0.707

#[2,] 0.707 -0.707

#

# Fraction of Variance 

#[1] 1 1

#

merge.forecast <- function (fx,n=1) {

x <- splice(fx$pred,fx$forecast[[n]])

colnames(x) <- seriesNames(fx$data$output)

return(x)

}

AIC <- function(model) {informationTestsCalculations(model)[3]}

require(dse)

require(matlab)

f <- matrix( c( 1.496523e-12,    0,  -1.325921e-14,

                       3.591688e-13,    0,   2.642115e-15,

                        0.000000000,  0.00000000, 1.0000000000)

                        ,byrow=TRUE,nrow=3,ncol=3)

g <- matrix( c(0.806718608,  2.549691,  0.6327740,

                       0.009381695, -1.547697, -0.2257979,

                       0.000000000,  0.000000,  0.0000000)

                        ,byrow=TRUE,nrow=3,ncol=3)

h <- eye(2,3)

k <- (f[,1:2,drop=FALSE])

NA_RE <- SS(F=f,H=h,K=k,G=g,z0=c(-2.8729652,  0.1631968,  1.0000000000),

              output.names=c("NA1","NA2"))

print(NA_RE)

is.SS(NA_RE)

stability(NA_RE)

tfplot(NA_RE <- simulate(NA_RE,sampleT=250,start=0,input=WRE.fx))

#NA_RE.data <- simulate(NA_RE,sampleT=20,noise=matrix(0,20,2))

m <- l(NA_RE,TSdata(output=NA_RE.data,input=WRE.data))

NA_RE.f <- forecast(m,horizon=500,conditioning.inputs=WRE.fx)

tfplot(NA_RE.f)

AIC(m)

shockDecomposition(toSSChol(m))