#

# Cut-and-Paste Code Below into Window Above to Run

# The MESOPOTAMIA MODEL (-10,000BCE - 0 BCE)

require(dse)

require(matlab)

#

# MEASUREMENT MATRIX (T-Q-N) (Growth) (Q-N)

#          N      Q      T

#[1,] -0.685 -0.692 0.2288

#[2,]  0.194  0.130 0.9724

#[3,] -0.702  0.710 0.0451

#

# Fraction of Variance 

#[1] 0.674 0.990 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]}

#

f <- matrix( c(1.044092872, -0.069114853, 0.07030666, -0.09916907,

                      0.019853110,  0.835132120, 0.04749995,  0.06906795,

                     -0.001702441, -0.009133031, 0.73592826, -0.00236148,

                      0.000000000,  0.000000000, 0.00000000,  1.00000000

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

#

# To Stabilize, Uncomment Next Line

# f[1,1] <- 0.9

#

h <- eye(3,4)

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

MESO <- SS(F=f,H=h,K=k,z0=c(-0.09916907,  0.06906795, -0.00236148,  1.0000000),

              output.names=c("MESO1","MESO2","MESO3"))

stability(MESO)

shockDecomposition(toSSChol(MESO))

#tfplot(simulate(MESO,sampleT=50,noise=matrix(0,50,2),start=1))

MESO.data <- simulate(MESO,sampleT=50,start=1)

m <- l(MESO,MESO.data)

#tfplot(m)

MESO.f <- forecast(m,horizon=50)

tfplot(MESO.f)

AIC(m)