diff --git a/R/dsemRTMB.R b/R/dsemRTMB.R
new file mode 100644
index 0000000..9c410e1
--- /dev/null
+++ b/R/dsemRTMB.R
@@ -0,0 +1,256 @@
+
+
+dsemRTMB <-
+function( sem,
+          tsdata,
+          family = rep("fixed",ncol(tsdata)),
+          estimate_delta0 = FALSE,
+          prior_negloglike = NULL,
+          control = dsem_control(),
+          covs = colnames(tsdata) ){
+
+  # General error checks
+  if( isFALSE(is(control, "dsem_control")) ) stop("`control` must be made by `dsem_control()`")
+  if( control$gmrf_parameterization=="projection" ){
+    if( any(family=="fixed" & colSums(!is.na(tsdata))>0) ){
+      stop("`family` cannot be `fixed` using `gmrf_parameterization=projection` for any variable with data")
+    }
+  }
+  if( isFALSE(is(tsdata,"ts")) ) stop("`tsdata` must be a `ts` object")
+
+  # General warnings
+  if( isFALSE(control$quiet) ){
+    tsdata_SD = apply( tsdata, MARGIN=2, FUN=sd, na.rm=TRUE )
+    if( any((max(tsdata_SD)/min(tsdata_SD)) > 10, rm.na=TRUE) ){
+      warning("Some variables in `tsdata` have much higher variance than others. Please consider rescaling variables to prevent issues with numerical convergence.")
+    }
+  }
+
+  # Build RAM
+  model = read_model( sem = sem,
+            times = as.numeric(time(tsdata)),
+            variables = colnames(tsdata),
+            covs = covs,
+            quiet = FALSE )
+  ram = make_matrices(
+            model = model,
+            times = as.numeric(time(tsdata)),
+            variables = colnames(tsdata) )
+
+  #
+  options = c(
+    ifelse(control$gmrf_parameterization=="separable", 0, 1),
+    switch(control$constant_variance, "conditional"=0, "marginal"=1, "diagonal"=2)
+  )
+  y_tj = tsdata
+
+  # Construct parameters
+  if( is.null(control$parameters) ){
+    Params = list(
+      #beta_p2 = subset(model,direction==2)$start,
+      #beta_p1 = subset(model,direction==1)$start,
+      beta_z = 0.01 * rnorm(n_p),
+      lnsigma_j = rep(0, n_j),
+      mu_j = rep(0, n_j),
+      delta0_j = rep(0, n_j),
+      x_tj = ifelse( is.na(tsdata), 0, tsdata )
+    )
+    #if( control$gmrf_parameterization=="separable" ){
+    #  Params$x_tj = ifelse( is.na(tsdata), 0, tsdata )
+    #}else{
+    #  Params$eps_tj = ifelse( is.na(tsdata), 0, tsdata )
+    #}
+
+    # Turn off initial conditions
+    if( estimate_delta0==FALSE ){
+      Params$delta0_j = numeric(0)
+    }
+
+    # Scale starting values with higher value for two-headed than one-headed arrows
+    which_nonzero = which(model[,5]>0)
+    beta_type = tapply( model[which_nonzero,8], INDEX=model[which_nonzero,5], max)
+    Params$beta_z = ifelse(beta_type==1, 0.01, 1)
+
+    # Override starting values if supplied
+    which_nonzero = which(model[,5]>0)
+    start_z = tapply( as.numeric(model[which_nonzero,4]), INDEX=model[which_nonzero,5], mean )
+    Params$beta_z = ifelse( is.na(start_z), Params$beta_z, start_z)
+  }else{
+    Params = control$parameters
+  }
+
+  # Construct map
+  if( is.null(control$map) ){
+    Map = list()
+    Map$x_tj = factor(ifelse( is.na(as.vector(tsdata)) | (family[col(tsdata)] %in% c("normal","binomial","poisson","gamma")), seq_len(prod(dim(tsdata))), NA ))
+    Map$lnsigma_j = factor( ifelse(family=="fixed", NA, seq_along(Params$lnsigma_j)) )
+
+    # Map off mean for latent variables
+    Map$mu_j = factor( ifelse(colSums(!is.na(tsdata))==0, NA, 1:ncol(tsdata)) )
+  }else{
+    Map = control$map
+  }
+
+  # Define random
+  if(isTRUE(control$use_REML)){
+    Random = c( "x_tj", "mu_j" )
+  }else{
+    Random = "x_tj"
+  }
+
+  # Build object
+  #source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "make_matrices.R") )
+  #source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "get_jnll.R") )
+  cmb <- function(f, ...) function(p) f(p, ...) ## Helper to make closure
+  #f(parlist, model, tsdata, family)
+  obj = RTMB::MakeADFun(
+          func = cmb( get_jnll,
+                      model = model,
+                      tsdata = tsdata,
+                      family = family ),
+          parameters = Params,
+          random = Random,
+          map = Map,
+          profile = control$profile,
+          silent = TRUE )
+
+  #
+  if( FALSE ){
+    repRTMB = obj$report()
+    range(Rep$Gamma_kk - repRTMB$Gamma_kk)
+    range(Rep$Rho_kk - repRTMB$Rho_kk)
+    range(Rep$Q_kk - repRTMB$Q_kk)
+    range(Rep$mu_tj - repRTMB$mu_tj)
+    range(Rep$delta_tj - repRTMB$delta_tj)
+    range(Rep$xhat_tj - repRTMB$xhat_tj)
+    range(Rep$jnll_gmrf - repRTMB$jnll_gmrf)
+
+    #
+    dgmrf( repRTMB$z_tj, mu=repRTMB$xhat_tj + repRTMB$delta_tj, Q=repRTMB$Q_kk, log=TRUE )
+    dgmrf( Rep$z_tj, mu=Rep$xhat_tj + Rep$delta_tj, Q=Rep$Q_kk, log=TRUE )
+    dgmrf( as.vector(Rep$z_tj), mu=as.vector(Rep$xhat_tj + Rep$delta_tj), Q=Rep$Q_kk, log=TRUE )
+    mvtnorm::dmvnorm( as.vector(Rep$z_tj - (Rep$xhat_tj + Rep$delta_tj)), sigma=as.matrix(solve(Rep$Q_kk)), log=TRUE )
+    dGMRF( as.vector(Rep$z_tj), mu=as.vector(Rep$xhat_tj + Rep$delta_tj), Q=Rep$Q_kk )
+
+    #
+    Rep = fit$obj$report( fit$obj$env$last.par )
+
+    #
+    example = list(
+      y = repRTMB$z_tj,
+      mu = repRTMB$xhat_tj + repRTMB$delta_tj,
+      Q = Rep$Q_kk
+    )
+    saveRDS( example, file=file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\scratch)','example.RDS') )
+  }
+
+  if(control$quiet==FALSE) list_parameters(obj)
+  # bundle
+  internal = list(
+    sem = sem,
+    tsdata = tsdata,
+    family = family,
+    estimate_delta0 = estimate_delta0,
+    control = control,
+    covs = covs
+  )
+
+  # Parse priors
+#  if( !is.null(prior_negloglike) ){
+#    # prior_negloglike = \(obj) -dnorm(obj$par[1],0,1,log=TRUE)
+#    prior_value = tryCatch( expr = prior_negloglike(obj) )
+#    if( is.na(prior_value) ) stop("Check `prior_negloglike(obj$par)`")
+#    obj$fn_orig = obj$fn
+#    obj$gr_orig = obj$gr
+#
+#    # BUild prior evaluator
+#    requireNamespace(RTMB)
+#    priors_obj = RTMB::MakeADFun( func = prior_negloglike,
+#                                  parameters = list(par=obj$par),
+#                                  silent = TRUE )
+#    obj$fn = \(pars) obj$fn_orig(pars) + priors_obj$fn(pars)
+#    obj$gr = \(pars) obj$gr_orig(pars) + priors_obj$gr(pars)
+#    internal$priors_obj = priors_obj
+#  }
+
+  # Further bundle
+  out = list( "obj" = obj,
+              "ram" = ram,
+              "sem_full" = model,
+              "tmb_inputs"=list("parameters"=Params, "random"=Random, "map"=Map),
+              #"call" = match.call(),
+              "internal" = internal )
+
+  # Export stuff
+  if( control$run_model==FALSE ){
+    return( out )
+  }
+
+  # Fit
+  #out$opt = fit_tmb( obj,
+  #                   quiet = control$quiet,
+  #                   control = list(eval.max=10000, iter.max=10000, trace=ifelse(control$quiet==TRUE,0,1) ),
+  #                   ... )
+
+  # Optimize
+  out$opt = list( "par"=obj$par )
+  for( i in seq_len(max(0,control$nlminb_loops)) ){
+    if( isFALSE(control$quiet) ) message("Running nlminb_loop #", i)
+    out$opt = nlminb( start = out$opt$par,
+                  objective = obj$fn,
+                  gradient = obj$gr,
+                  upper = control$upper,
+                  lower = control$lower,
+                  control = list( eval.max = control$eval.max,
+                                  iter.max = control$iter.max,
+                                  trace = control$trace ) )
+  }
+
+  # Newtonsteps
+  for( i in seq_len(max(0,control$newton_loops)) ){
+    if( isFALSE(control$quiet) ) message("Running newton_loop #", i)
+    g = as.numeric( obj$gr(out$opt$par) )
+    h = optimHess(out$opt$par, fn=obj$fn, gr=obj$gr)
+    out$opt$par = out$opt$par - solve(h, g)
+    out$opt$objective = obj$fn(out$opt$par)
+  }
+  out$internal$parhat = obj$env$parList()
+
+  if( isTRUE(control$extra_convergence_checks) ){
+    # Gradient checks
+    Grad_fixed = obj$gr( out$opt$par )
+    if( isTRUE(any(Grad_fixed > 0.01)) ){
+      warning("Some gradients are higher than 0.01. Some parameters might not be converged.  Consider increasing `control$newton_loops`")
+    }
+    # Hessian check ... condition and positive definite
+    Hess_fixed = optimHess( par=out$opt$par, fn=obj$fn, gr=obj$gr, control=list(ndeps=rep(0.001,length(out$opt$par))) )
+    Eigen_fixed = eigen( Hess_fixed, only.values=TRUE )
+    if( (max(Eigen_fixed$values)/min(Eigen_fixed$values)) > 1e6 ){
+      # See McCullough and Vinod 2003
+      warning("The ratio of maximum and minimum Hessian eigenvalues is high. Some parameters might not be identifiable.")
+    }
+    if( isTRUE(any(Eigen_fixed$values < 0)) ){
+      warning("Some Hessian eigenvalue is negative. Some parameters might not be converged.")
+    }
+  }else{
+    Hess_fixed = NULL
+  }
+
+  # Run sdreport
+  if( isTRUE(control$getsd) ){
+    if( isTRUE(control$verbose) ) message("Running sdreport")
+    if( is.null(Hess_fixed) ){
+      Hess_fixed = optimHess( par=out$opt$par, fn=obj$fn, gr=obj$gr, control=list(ndeps=rep(0.001,length(out$opt$par)))  )
+    }
+    out$sdrep = TMB::sdreport( obj,
+                          par.fixed = out$opt$par,
+                          hessian.fixed = Hess_fixed,
+                          getJointPrecision = control$getJointPrecision )
+  }else{
+    out$sdrep = NULL
+  }
+
+  # output
+  class(out) = "dsem"
+  return(out)
+}
diff --git a/R/get_jnll.R b/R/get_jnll.R
new file mode 100644
index 0000000..72d1b86
--- /dev/null
+++ b/R/get_jnll.R
@@ -0,0 +1,153 @@
+# model function
+get_jnll <-
+function( parlist,
+          model,
+          tsdata,
+          family ) {
+
+  "c" <- ADoverload("c")
+  "[<-" <- ADoverload("[<-")
+
+  # Unpack parameters explicitly
+  beta_z = parlist$beta_z
+  delta0_j = parlist$delta0_j
+  mu_j = parlist$mu_j
+  sigma_j = exp( parlist$lnsigma_j )
+  x_tj = parlist$x_tj
+
+  #
+  n_z = length(unique(model$parameter))
+  #n_p2 = length(unique(subset(model,direction==2)$parameter))
+  #n_p1 = length(unique(subset(model,direction==1)$parameter))
+  n_t = nrow(tsdata)
+  n_j = ncol(tsdata)
+  n_k = prod(dim(tsdata))
+
+  # Unpack
+  #model_unique = model[match(unique(model$parameter),model$parameter),]
+  #beta_p[which(model_unique$direction==1)] = beta_p1
+  #beta_p[which(model_unique$direction==2)] = beta_p2
+
+  # Build matrices
+  ram = make_matrices(
+            beta_p = beta_z,
+            model = model,
+            times = as.numeric(time(tsdata)),
+            variables = colnames(tsdata) )
+
+  # Assemble
+  Rho_kk = AD(ram$P_kk)
+  IminusRho_kk = Diagonal(n_k) - Rho_kk
+  # Assemble variance
+  Gamma_kk = invV_kk = AD(ram$G_kk)
+  invV_kk@x = 1 / Gamma_kk@x^2
+
+  # Calculate effect of initial condition -- SPARSE version
+  delta_tj = matrix( 0, nrow=n_t, ncol=n_j )
+  if( length(delta0_j)>0 ){
+    delta_tj[1,] = delta0_j
+    delta_k = solve( IminusRho_kk, as.vector(delta_tj) )
+    delta_tj[] = delta_k
+  }
+
+  #
+  xhat_tj = outer( rep(1,n_t), mu_j )
+  z_tj = x_tj
+
+  # Probability of GMRF
+  Q_kk = t(IminusRho_kk) %*% invV_kk %*% IminusRho_kk
+  jnll_gmrf = -1 * dgmrf( as.vector(z_tj), mu=as.vector(xhat_tj + delta_tj), Q=Q_kk, log=TRUE )
+  REPORT( Q_kk )
+
+  # Likelihood of data | random effects
+  loglik_tj = mu_tj = devresid_tj = matrix( 0, nrow=n_t, ncol=n_j )
+  pow = function(a,b) a^b
+  for( t in 1:n_t ){
+  for( j in 1:n_j ){
+    # familycode = 0 :  don't include likelihood
+    if( family[j]=="fixed" ){
+      mu_tj[t,j] = z_tj[t,j];
+      if(!is.na(y_tj[t,j])){
+        #SIMULATE{
+        #  y_tj[t,j] = mu_tj[t,j];
+        #}
+      }
+      devresid_tj[t,j] = 0;
+    }
+    # familycode = 1 :  normal
+    if( family[j]=="normal" ){
+      mu_tj[t,j] = z_tj[t,j];
+      if(!is.na(y_tj[t,j])){
+        loglik_tj[t,j] = dnorm( y_tj[t,j], mu_tj[t,j], sigma_j[j], TRUE );
+        #SIMULATE{
+        #  y_tj[t,j] = rnorm( mu_tj[t,j], sigma_j[j] );
+        #}
+      }
+      devresid_tj[t,j] = y_tj[t,j] - mu_tj[t,j];
+    }
+    # familycode = 2 :  binomial
+    if( family[j]=="binomial" ){
+      mu_tj[t,j] = invlogit(z_tj[t,j]);
+      if(!is.na(y_tj[t,j])){
+        loglik_tj[t,j] = dbinom( y_tj[t,j], Type(1.0), mu_tj[t,j], TRUE );
+        #SIMULATE{
+        #  y_tj[t,j] = rbinom( Type(1), mu_tj[t,j] );
+        #}
+      }
+      devresid_tj[t,j] = sign(y_tj[t,j] - mu_tj[t,j]) * pow(-2*(((1-y_tj[t,j])*log(1-mu_tj[t,j]) + y_tj[t,j]*log(mu_tj[t,j]))), 0.5);
+    }
+    # familycode = 3 :  Poisson
+    if( family[j]=="poisson" ){
+      mu_tj[t,j] = exp(z_tj[t,j]);
+      if(!is.na(y_tj[t,j])){
+        loglik_tj[t,j] = dpois( y_tj[t,j], mu_tj[t,j], TRUE );
+        #SIMULATE{
+        #  y_tj[t,j] = rpois( mu_tj[t,j] );
+        #}
+      }
+      devresid_tj[t,j] = sign(y_tj[t,j] - mu_tj[t,j]) * pow(2*(y_tj[t,j]*log((Type(1e-10) + y_tj[t,j])/mu_tj[t,j]) - (y_tj[t,j]-mu_tj[t,j])), 0.5);
+    }
+    # familycode = 4 :  Gamma:   shape = 1/CV^2; scale = mean*CV^2
+    if( family[j]=="gamma" ){
+      mu_tj[t,j] = exp(z_tj[t,j]);
+      if(!is.na(y_tj[t,j])){
+        loglik_tj[t,j] = dgamma( y_tj[t,j], pow(sigma_j[j],-2), mu_tj[t,j]*pow(sigma_j[j],2), TRUE );
+        #SIMULATE{
+        #  y_tj[t,j] = rgamma( pow(sigma_j[j],-2), mu_tj[t,j]*pow(sigma_j[j],2) );
+        #}
+      }
+      devresid_tj[t,j] = sign(y_tj[t,j] - mu_tj[t,j]) * pow(2 * ( (y_tj[t,j]-mu_tj[t,j])/mu_tj[t,j] - log(y_tj[t,j]/mu_tj[t,j]) ), 0.5);
+    }
+  }}
+  jnll = -1 * sum(loglik_tj);
+  jnll = jnll + jnll_gmrf;
+
+  #
+  REPORT( loglik_tj )
+  REPORT( jnll_gmrf )
+  REPORT( xhat_tj ); # needed to simulate new GMRF in R
+  #REPORT( delta_k ); # FIXME>  Eliminate in simulate.dsem
+  REPORT( delta_tj ); # needed to simulate new GMRF in R
+  REPORT( Rho_kk );
+  REPORT( Gamma_kk );
+  REPORT( mu_tj );
+  REPORT( devresid_tj );
+  REPORT( IminusRho_kk );
+  REPORT( jnll );
+  REPORT( loglik_tj );
+  REPORT( jnll_gmrf );
+  #SIMULATE{
+  #  REPORT( y_tj );
+  #}
+  REPORT( z_tj );
+  ADREPORT( z_tj );
+
+  jnll
+}
+
+if( FALSE ){
+  family = rep("fixed",ncol(tsdata))
+  y_tj = tsdata
+  estimate_delta0 = FALSE
+}
+
diff --git a/R/make_matrices.R b/R/make_matrices.R
new file mode 100644
index 0000000..05ff3b1
--- /dev/null
+++ b/R/make_matrices.R
@@ -0,0 +1,55 @@
+make_matrices <-
+function( beta_p,
+          model,
+          times,
+          variables ){
+
+  if(missing(beta_p)){
+    model_unique = model[match(unique(model$parameter),model$parameter),]
+    beta_p =  as.numeric(model_unique$start)
+  }
+
+  # Loop through paths
+  P_kk = drop0(sparseMatrix( i=1, j=1, x=0, dims=rep(length(variables)*length(times),2) ))   # Make with a zero
+  #P_kk = AD(P_kk)
+  G_kk = (P_kk)
+  for( i in seq_len(nrow(model)) ){
+    lag = as.numeric(model[i,2])
+    L_tt = sparseMatrix( i = seq(lag+1,length(times)),
+                         j = seq(1,length(times)-lag),
+                         x = 1,
+                         dims = rep(length(times),2) )
+
+    P_jj = sparseMatrix( i = match(model[i,'second'],variables),
+                         j = match(model[i,'first'],variables),
+                         x = 1,
+                         dims = rep(length(variables),2) )
+
+    # Assemble
+    if(abs(as.numeric(model[i,'direction']))==1){
+      tmp_kk = (kronecker(P_jj, L_tt))
+      P_kk = P_kk + beta_p[model$parameter[i]] * tmp_kk # AD(tmp_kk)
+    }else{
+      tmp_kk = (kronecker(P_jj, L_tt))
+      G_kk = G_kk + beta_p[model$parameter[i]] * tmp_kk # AD(tmp_kk)
+    }
+  }
+
+  # Diagonal component
+  I_kk = Diagonal(nrow(P_kk))
+
+  # Assemble
+  IminusP_kk = I_kk - P_kk
+  invV_kk = AD(G_kk)
+  invV_kk@x = 1 / G_kk@x^2
+  Q_kk = t(IminusP_kk) %*% invV_kk %*% IminusP_kk
+
+  out = list(
+    "P_kk" = P_kk,
+    "G_kk" = G_kk,
+    "invV_kk" = invV_kk,
+    "IminusP_kk" = IminusP_kk,
+    "Q_kk" = Q_kk
+  )
+  return(out)
+}
diff --git a/R/read_model.R b/R/read_model.R
new file mode 100644
index 0000000..c0a44b2
--- /dev/null
+++ b/R/read_model.R
@@ -0,0 +1,138 @@
+
+#' @title Make a RAM (Reticular Action Model)
+#'
+#' @description \code{read_model} converts SEM arrow notation to \code{model} describing SEM parameters
+#'
+#' @inheritParams dsem
+#' @param times A character vector listing the set of times in order
+#' @param variables A character vector listing the set of variables
+#' @param covs A character vector listing variables for which to estimate a standard deviation
+#' @param quiet Boolean indicating whether to print messages to terminal
+#'
+#' @details
+#' \strong{RAM specification using arrow-and-lag notation}
+#'
+#'
+#' @export
+read_model <-
+function( sem,
+          times,
+          variables,
+          covs = NULL,
+          quiet = FALSE ){
+
+  "c" <- ADoverload("c")
+  "[<-" <- ADoverload("[<-")
+
+  ####### Error checks
+  if( !is.numeric(times) ) stop("`times` must be numeric in `make_dsem_ram`")
+
+  ####### Define local functions
+  #
+  add.variances <- function() {
+      variables <- need.variance()
+      nvars <- length(variables)
+      if (nvars == 0)
+          return(model)
+      message("NOTE: adding ", nvars, " variances to the model")
+      paths <- character(nvars)
+      par.names <- character(nvars)
+      for (i in 1:nvars) {
+          paths[i] <- paste(variables[i], "<->", variables[i])
+          par.names[i] <- paste("V[", variables[i], "]", sep = "")
+      }
+      model.2 <- cbind(
+        'path' = c(model[, 1], paths),
+        'lag' = c(model[,2], rep(0,nvars)),
+        'name' = c(model[, 3], par.names),
+        'start' = c(model[, 4], rep(NA, length(paths))) )
+      model.2
+  }
+  need.variance <- function() {
+      all.vars <- classify_variables(model)
+      exo.vars <- all.vars$exogenous
+      end.vars <- all.vars$endogenous
+      variables <- logical(0)
+      for (i in seq_len(nrow(model))) {
+          paths = model[i,1]
+          lag = model[i,2]
+          vars <- gsub(pattern=" ", replacement="", x=paths)
+          vars <- sub("-*>", "->", sub("<-*", "<-", vars))
+          vars <- sub("<->|<-", "->", vars)
+          vars <- strsplit(vars, "->")[[1]]
+          if ((vars[1] != vars[2]) | (lag != 0)) {
+              for (a.variable in vars) {
+                if (is.na(variables[a.variable]))
+                  variables[a.variable] <- TRUE
+              }
+          }
+          else {
+              variables[vars[1]] <- FALSE
+          }
+      }
+      if (!exog.variances && length(exo.vars) > 0)
+          variables[exo.vars] <- FALSE
+      if (!endog.variances && length(end.vars) > 0)
+          variables[end.vars] <- FALSE
+      names(variables)[variables]
+  }
+
+  ####### Step 2 -- Make RAM
+  # convert to data frame
+  model = scan( text = sem,
+                what = list(path = "", lag = 1, par = "", start = 1, dump = ""),
+                sep = ",",
+                strip.white = TRUE,
+                comment.char = "#",
+                fill = TRUE,
+                quiet = quiet)
+  model$path <- gsub("\\t", " ", model$path)
+  model$par[model$par == ""] <- NA
+  model <- data.frame( "path"=model$path, "lag"=model$lag,
+                       "name"=model$par, "start"=model$start)
+
+  # Adding a SD automatically
+  if( !is.null(covs) ){
+    for (cov in covs) {
+      vars <- strsplit(cov, "[ ,]+")[[1]]
+      nvar <- length(vars)
+      for (i in 1:nvar) {
+      for (j in i:nvar) {
+        p1 = paste(vars[i], "<->", vars[j])
+        p2 = if (i==j) paste("V[", vars[i], "]", sep = "") else paste("C[",vars[i], ",", vars[j], "]", sep = "")
+        p3 = NA
+        row <- c(p1, 0, p2, p3)
+        if( any((row[1]==model[,1]) & (row[2]==model[,2])) ){
+          next
+        }else{
+          model <- rbind(model, row, deparse.level = 0)
+        }
+      }}
+    }
+  }
+
+  exog.variances = endog.variances = TRUE
+  model = add.variances()
+
+  # Add parameter column
+  par.names = model[, 3]
+  pars = na.omit(unique(par.names))
+  par.nos = apply(outer(pars, par.names, "=="), 2, which)
+  par.nos = unlist(sapply( par.nos, FUN=\(x) ifelse(length(x)==0, 0, x) ))
+  model = cbind( model, "parameter"=par.nos )
+
+  # Add incidence to model
+  model = cbind( model, first=NA, second=NA, direction=NA )
+  for( i in seq_len(nrow(model)) ){
+    path = parse_path(model[i,1])
+    model[i,c('first','second','direction')] = unlist( path[c('first','second','direction')] )
+  }
+
+  ####### Step 2 -- Make RAM
+
+  # Deal with fixed values
+  #if( max(model$parameter) != length(beta_p) ) stop("Check beta_p")
+
+  return(model)
+}
+
diff --git a/scratch/example.R b/scratch/example.R
new file mode 100644
index 0000000..70b3b49
--- /dev/null
+++ b/scratch/example.R
@@ -0,0 +1,24 @@
+
+
+setwd( R'(C:\Users\James.Thorson\Desktop\Git\dsem\scratch)' )
+
+example = readRDS( "example.RDS" )
+attach(example)
+
+# Custom version
+dGMRF = function(x, mu, Q ){
+  #(2*pi)^(-n_k/2) *
+  #determinant(solve(Q))^(-1/2) *
+  #exp(-1/2 * (x-mu) %*% Q %*% (x-mu) )
+  (-length(as.vector(x))/2) * log(2*pi) +
+  (-1/2) * -1 * determinant(Q, log=TRUE)$modulus +
+  (-1/2 * as.vector(x-mu) %*% Q %*% as.vector(x-mu) )[1,1]
+}
+
+#
+library(RTMB)
+dgmrf( y, mu=mu, Q=Q, log=TRUE )
+dgmrf( as.vector(y), mu=as.vector(mu), Q=Q, log=TRUE )
+mvtnorm::dmvnorm( y, mu, sigma=as.matrix(solve(Q)), log=TRUE )
+mvtnorm::dmvnorm( as.vector(y), as.vector(mu), sigma=as.matrix(solve(Q)), log=TRUE )
+dGMRF( y, mu=mu, Q=Q )
diff --git a/scratch/example.RDS b/scratch/example.RDS
new file mode 100644
index 0000000..56452ce
Binary files /dev/null and b/scratch/example.RDS differ
diff --git a/scratch/test_dsemRTMB.R b/scratch/test_dsemRTMB.R
new file mode 100644
index 0000000..6b2e08c
--- /dev/null
+++ b/scratch/test_dsemRTMB.R
@@ -0,0 +1,49 @@
+
+library(dsem)
+library(RTMB)
+
+# Define model
+sem = "
+  # Link, lag, param_name
+  cprofits -> consumption, 0, a1
+  cprofits -> consumption, 1, a2
+  pwage -> consumption, 0, a3
+  gwage -> consumption, 0, a3
+  cprofits -> invest, 0, b1
+  cprofits -> invest, 1, b2
+  capital -> invest, 0, b3
+  gnp -> pwage, 0, c2
+  gnp -> pwage, 1, c3
+  time -> pwage, 0, c1
+"
+
+# Load data
+data(KleinI, package="AER")
+TS = ts(data.frame(KleinI, "time"=time(KleinI) - 1931))
+tsdata = TS[,c("time","gnp","pwage","cprofits",'consumption',
+               "gwage","invest","capital")]
+
+# Fit model
+fit = dsem( sem=sem,
+            tsdata = tsdata,
+            estimate_delta0 = TRUE,
+            control = dsem_control( quiet=TRUE ) )
+#ParHat = fit$internal$parhat
+#Rep = fit$obj$report()
+
+# RUN dsemRTMB line-by-line
+if( FALSE ){
+  control = dsem_control()
+  covs = colnames(tsdata)
+}
+
+#
+source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "make_matrices.R") )
+source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "get_jnll.R") )
+source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "read_model.R") )
+source( file.path(R'(C:\Users\James.Thorson\Desktop\Git\dsem\R)', "dsemRTMB.R") )
+fitRTMB = dsemRTMB( sem = sem,
+            tsdata = tsdata,
+            estimate_delta0 = TRUE,
+            control = dsem_control( quiet = TRUE) )
+