Changes made in revision 4 of KalmanEM.r

--- /tmp/cvx_GdtGgM	2012-02-07 19:42:31.000000000 -0800
+++ /tmp/cvx_lCjzBv	2012-02-07 19:42:31.000000000 -0800
@@ -1,7 +1,7 @@
 KalmanEM = function(y, whichPop=NA, Binit=NA, Uinit=0, Qinit=0.05, Ainit=0, Rinit=0.05, x00=NA, V00init=10, Zinit=NA, max.iter=5000, varcov.Q="unconstrained", varcov.R="diagonal", B.constraint="identity", U.constraint="unconstrained", Z.constraint="missingvalues", A.groups=NA, U.groups=NA, Q.groups=NA, R.groups=NA, x00.groups=NA, B.groups=NA, U.bounds=c(-1,1), logQ.bounds = c(log(1.0e-05),log(1.0)), logR.bounds = c(log(1.0e-05),log(1.0)), x00prior=FALSE, estInteractions=FALSE, MonteCarloInit = FALSE, numInits = 500, numInitSteps = 10, tol=0.01, silent=FALSE, debugit=FALSE) 
 {
 #Version
-vers = "2.67"  #EEH July 17, 2009
+vers = "2.68"  #EEH Sept 24, 2009
                                                                                                                                     #comment on x00prior=FALSE; x00prior=TRUE case is not programmed in yet.  In that case x00 fixed, not estimated.
 #x00prior=FALSE means that we treat x00 as an unknown fixed parameter that is estimated.  So P_0=V00=0, but we don't want to set
 #V00=0 because then the EM algorithm won't converge.  We want V00 big, and then at the end, return the likelihood with V00=0                                                                                                                                          
@@ -467,7 +467,7 @@
             }
          else {
             stop("Currently Kalman-EM can't estimate B diagonal if Q not diagonal too")
-            # EEH: I have figured out the update equation for this case.
+            # EEH: I have not figured out the update equation for this case.
             }
          
          }
@@ -1284,21 +1284,49 @@
 TIC = try(-2*mssm.model$loglike+2*sum(diag(J0%*%solve(I0))))
 mssm.model$TIC = ifelse(!inherits(TIC, "try-error"),TIC,NA)
 
-#standard errors
-stderr=try(sqrt(diag(solve(emhess$Hessian))))      #invert the Hessian; wrap in a try() in case it fails
-stderr=ifelse(!inherits(stderr, "try-error"),stderr,rep(NA,length(paramvector)))
-stderr.model=vectorizemssm(mssm.model,paramvector=stderr)
-if(mssm.model$x00prior==TRUE) {stderr.model$x00=NA; stderr.model$V00=NA}
-if(mssm.model$x00prior==FALSE) stderr.model$V00=NA
-if(mssm.model$estInteractions==FALSE) stderr.model$B=NA
-if(mssm.model$estInteractions==TRUE) stderr.model$B="Hessian when B estimated is not coded in yet"
-mssm.model$A.stderr=stderr.model$A
-mssm.model$B.stderr=stderr.model$B
-mssm.model$U.stderr=stderr.model$U
-mssm.model$Q.stderr=stderr.model$Q
-mssm.model$R.stderr=stderr.model$R
-mssm.model$x00.stderr=stderr.model$x00
-mssm.model$V00.stderr=stderr.model$V00
+return(mssm.model)
+}
 
+#######################################################################################################
+#   computeCIs
+#######################################################################################################
+computeCIs = function(mssm.model, method="Hessian", alpha=0.05) {
+if(method=="Hessian")  {
+    #if the model has no Hessian specified, then run emHessian to get it
+    if(is.null(mssm.model$Hessian)) mssm.model=emHessian(mssm.model)
+    #standard errors
+    stderr=try(sqrt(diag(solve(mssm.model$Hessian))))      #invert the Hessian; wrap in a try() in case it fails
+    if(inherits(stderr, "try-error")) {
+        stderr=rep(NA,length(paramvector))
+        warning("Hessian cannot be inverted; stderr = NA is being returned")
+        }
+    stderr.model=vectorizemssm(mssm.model,paramvector=stderr)
+    if(mssm.model$x00prior==TRUE) {stderr.model$x00=NA; stderr.model$V00=NA}
+    if(mssm.model$x00prior==FALSE) stderr.model$V00=NA
+    if(mssm.model$estInteractions==FALSE) stderr.model$B=NA
+    if(mssm.model$estInteractions==TRUE) stderr.model$B="Hessian when B estimated is not coded in yet"
+    mssm.model$A.stderr=stderr.model$A
+    mssm.model$B.stderr=stderr.model$B
+    mssm.model$U.stderr=stderr.model$U
+    mssm.model$Q.stderr=stderr.model$Q
+    mssm.model$R.stderr=stderr.model$R
+    mssm.model$x00.stderr=stderr.model$x00
+    mssm.model$V00.stderr=stderr.model$V00
+    mssm.model$A.CI.upper=mssm.model$A + qnorm(1-alpha/2)*stderr.model$A
+    mssm.model$A.CI.lower=mssm.model$A - qnorm(1-alpha/2)*stderr.model$A
+    mssm.model$B.CI.upper=mssm.model$A + qnorm(1-alpha/2)*stderr.model$B
+    mssm.model$B.CI.lower=mssm.model$A - qnorm(1-alpha/2)*stderr.model$B
+    mssm.model$U.CI.upper=mssm.model$U + qnorm(1-alpha/2)*stderr.model$U
+    mssm.model$U.CI.lower=mssm.model$U - qnorm(1-alpha/2)*stderr.model$U
+    mssm.model$Q.CI.upper=mssm.model$Q + qnorm(1-alpha/2)*stderr.model$Q
+    mssm.model$Q.CI.lower=mssm.model$Q - qnorm(1-alpha/2)*stderr.model$Q
+    mssm.model$R.CI.upper=mssm.model$R + qnorm(1-alpha/2)*stderr.model$R
+    mssm.model$R.CI.lower=mssm.model$R - qnorm(1-alpha/2)*stderr.model$R
+    mssm.model$x00.CI.upper=mssm.model$x00 + qnorm(1-alpha/2)*stderr.model$x00
+    mssm.model$x00.CI.lower=mssm.model$x00 - qnorm(1-alpha/2)*stderr.model$x00
+    mssm.model$V00.CI.upper=mssm.model$V00 + qnorm(1-alpha/2)*stderr.model$V00
+    mssm.model$V00.CI.lower=mssm.model$V00 - qnorm(1-alpha/2)*stderr.model$V00
+}
+else stop("Only method Hessian programmed currently")
 return(mssm.model)
-}
\ No newline at end of file
+}