diff --git a/.travis.yml b/.travis.yml
index 3d1fe0b..8e71d26 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,4 +1,5 @@
 language: r
+cache: packages
 sudo: required
 warnings_are_errors: false
 
diff --git a/NAMESPACE b/NAMESPACE
index f44d800..c62871a 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,18 +1,20 @@
 useDynLib(BMR)
 
-importFrom(Rcpp, evalCpp)
 import(ggplot2)
+importFrom(Rcpp, evalCpp)
 importFrom(grid, grid.newpage, viewport, pushViewport, grid.layout)
 importFrom("grDevices", "cairo_ps", "dev.list", "dev.off")
-importFrom("stats", "acf", "dbeta", "dgamma", "dnorm", "dunif",
-           "embed", "optim", "pacf", "qnorm", "rnorm", "runif", "sd",
-           "var")
+importFrom("stats", "acf", "embed", "optim", "pacf", "qnorm", "rnorm", "runif", "sd", "var")
 
 exportClasses("Rcpp_bvarm")
 exportClasses("Rcpp_bvars")
 exportClasses("Rcpp_bvars")
 exportClasses("Rcpp_cvar")
 exportClasses("Rcpp_bvartvp")
+exportClasses("Rcpp_gensys")
+exportClasses("Rcpp_uhlig")
+exportClasses("Rcpp_dsge_gensys")
+exportClasses("Rcpp_dsgevar_gensys")
 
 export(IRF,forecast,modecheck,prior,stationarity,gtsplot,gacf,gpacf)
 
@@ -22,13 +24,13 @@ S3method(plot, Rcpp_bvarm)
 S3method(plot, Rcpp_bvars)
 S3method(plot, Rcpp_bvarw)
 S3method(plot, Rcpp_bvartvp)
-S3method(plot, EDSGE)
-S3method(plot, DSGEVAR)
+S3method(plot, Rcpp_dsge_gensys)
+S3method(plot, Rcpp_dsgevar_gensys)
 
 S3method(forecast, Rcpp_bvarm)
 S3method(forecast, Rcpp_bvars)
 S3method(forecast, Rcpp_bvarw)
 S3method(forecast, Rcpp_cvar)
 
-S3method(modecheck, EDSGE)
-S3method(modecheck, DSGEVAR)
+S3method(modecheck, Rcpp_dsge_gensys)
+S3method(modecheck, Rcpp_dsgevar_gensys)
diff --git a/R/Plot.R b/R/Plot.R
index 5524d38..5964a9f 100644
--- a/R/Plot.R
+++ b/R/Plot.R
@@ -592,6 +592,6 @@ plot.Rcpp_dsgevar_gensys <- function(x,parnames=NULL,BinDenom=40,MCMCplot=FALSE,
 
 .plotdsgevar <- function(obj,parnames=NULL,BinDenom=40,trace_plot=FALSE,save=FALSE,height=13,width=13)
 {    
-    dsge_obj <- obj$get_dsge()
+    dsge_obj <- obj$dsge
     .plotedsge(dsge_obj,parnames,BinDenom,trace_plot,save,height,width)
 }
diff --git a/src/bmlib/include/dsge/dsgevar_class.tpp b/src/bmlib/include/dsge/dsgevar_class.tpp
index 525f00c..da4d17b 100644
--- a/src/bmlib/include/dsge/dsgevar_class.tpp
+++ b/src/bmlib/include/dsge/dsgevar_class.tpp
@@ -458,7 +458,7 @@ dsgevar<T>::IRF(const int n_irf_periods)
         arma::mat dsge_obs_impact = dsge_obj_copy.kalman_mat_H.t() * G_state*arma::chol(shocks_cov,"lower");
 
         arma::mat Q_dcm, R_dcm;
-        arma::qr(Q_dcm,R_dcm,dsge_obs_impact);
+        arma::qr(Q_dcm,R_dcm,dsge_obs_impact.t());
 
         Q_dcm *= -1.0;
         R_dcm *= -1.0;
diff --git a/src/modules/dsge_R.hpp b/src/modules/dsge_R.hpp
index f5845ce..d075147 100644
--- a/src/modules/dsge_R.hpp
+++ b/src/modules/dsge_R.hpp
@@ -64,6 +64,9 @@ class dsge_gensys_R : public bm::dsge<bm::gensys>
 
         void model_fn_R(const arma::vec& pars_inp, bm::gensys& lrem_obj_inp, arma::mat& shocks_cov_out, arma::mat& C_out, arma::mat& H_out, arma::mat& R_out);
 
+        //
+
+        SEXP get_model_fn();
         void set_model_fn(SEXP model_fn_inp);
         void eval_model(Rcpp::NumericVector pars_inp);
 
@@ -93,14 +96,15 @@ class dsgevar_gensys_R : public bm::dsgevar<bm::gensys>
         arma::vec mcmc_initial_lb;
         arma::vec mcmc_initial_ub;
 
+        void model_fn_R(const arma::vec& pars_inp, bm::gensys& lrem_obj_inp, arma::mat& shocks_cov_out, arma::mat& C_out, arma::mat& H_out, arma::mat& R_out);
+
         // member functions
 
         void build_R(const arma::mat& data_raw, bool cons_term_inp, int p_inp, double lambda_inp);
 
         arma::mat get_dsge_draws();
 
-        void model_fn_R(const arma::vec& pars_inp, bm::gensys& lrem_obj_inp, arma::mat& shocks_cov_out, arma::mat& C_out, arma::mat& H_out, arma::mat& R_out);
-
+        SEXP get_model_fn();
         void set_model_fn(SEXP model_fn_inp);
         void eval_model(Rcpp::NumericVector pars_inp);
 
diff --git a/src/modules/dsge_gensys_R.cpp b/src/modules/dsge_gensys_R.cpp
index 4beca39..8b7e208 100644
--- a/src/modules/dsge_gensys_R.cpp
+++ b/src/modules/dsge_gensys_R.cpp
@@ -53,12 +53,13 @@ RCPP_MODULE(dsge_gensys_module)
         .field( "mcmc_initial_lb", &dsge_gensys_R::mcmc_initial_lb )
         .field( "mcmc_initial_ub", &dsge_gensys_R::mcmc_initial_ub )
 
+        .property( "lrem", &dsge_gensys_R::get_lrem_R, &dsge_gensys_R::set_lrem_R )
+        // .property( "model", &dsge_gensys_R::get_model_fn, &dsge_gensys_R::set_model_fn )
+
+        .method( "get_model_fn", &dsge_gensys_R::get_model_fn )
         .method( "set_model_fn", &dsge_gensys_R::set_model_fn )
         .method( "eval_model", &dsge_gensys_R::eval_model )
 
-        .method( "get_lrem", &dsge_gensys_R::get_lrem_R )
-        .method( "set_lrem", &dsge_gensys_R::set_lrem_R )
-
         .method( "set_bounds", &dsge_gensys_R::set_bounds_R )
         .method( "set_prior", &dsge_gensys_R::set_prior_R )
 
@@ -72,6 +73,12 @@ RCPP_MODULE(dsge_gensys_module)
 //
 // wrapper functions to catch errors and handle memory pointers
 
+SEXP
+dsge_gensys_R::get_model_fn()
+{
+    return model_fn_SEXP;
+}
+
 void
 dsge_gensys_R::set_model_fn(SEXP model_fn_inp)
 {
@@ -92,6 +99,8 @@ dsge_gensys_R::set_model_fn(SEXP model_fn_inp)
     }
 }
 
+//
+
 void
 dsge_gensys_R::eval_model(Rcpp::NumericVector pars_inp)
 {
@@ -151,14 +160,18 @@ dsge_gensys_R::model_fn_R(const arma::vec& pars_inp, bm::gensys& lrem_obj_inp, a
     }
 }
 
-gensys_R dsge_gensys_R::get_lrem_R()
+// get and set LREM
+
+gensys_R
+dsge_gensys_R::get_lrem_R()
 {
     gensys_R lrem_obj_out = static_cast<gensys_R&>(lrem_obj);
 
     return lrem_obj_out;
 }
 
-void dsge_gensys_R::set_lrem_R(gensys_R lrem_obj_inp)
+void
+dsge_gensys_R::set_lrem_R(gensys_R lrem_obj_inp)
 {
     try {
         lrem_obj = static_cast<bm::gensys&>(lrem_obj_inp);
@@ -171,7 +184,8 @@ void dsge_gensys_R::set_lrem_R(gensys_R lrem_obj_inp)
 
 // set bounds and prior
 
-void dsge_gensys_R::set_bounds_R(arma::vec lower_bounds_inp, arma::vec upper_bounds_inp)
+void
+dsge_gensys_R::set_bounds_R(arma::vec lower_bounds_inp, arma::vec upper_bounds_inp)
 {
     try {
         const int n_vals = lower_bounds_inp.n_elem;
@@ -194,7 +208,8 @@ void dsge_gensys_R::set_bounds_R(arma::vec lower_bounds_inp, arma::vec upper_bou
     }
 }
 
-void dsge_gensys_R::set_prior_R(const arma::uvec& prior_form_inp, const arma::mat& prior_pars_inp)
+void
+dsge_gensys_R::set_prior_R(const arma::uvec& prior_form_inp, const arma::mat& prior_pars_inp)
 { 
     try {
         this->set_prior(prior_form_inp,prior_pars_inp);
@@ -207,7 +222,8 @@ void dsge_gensys_R::set_prior_R(const arma::uvec& prior_form_inp, const arma::ma
 
 // 
 
-SEXP dsge_gensys_R::estim_mode_R(const arma::vec& initial_vals)
+SEXP
+dsge_gensys_R::estim_mode_R(const arma::vec& initial_vals)
 {
     try {
         optim::opt_settings settings;
@@ -231,7 +247,8 @@ SEXP dsge_gensys_R::estim_mode_R(const arma::vec& initial_vals)
     return R_NilValue;
 }
 
-void dsge_gensys_R::estim_mcmc_R(const arma::vec& initial_vals, int n_pop, int n_gen, int n_burnin)
+void
+dsge_gensys_R::estim_mcmc_R(const arma::vec& initial_vals, int n_pop, int n_gen, int n_burnin)
 {
     try {
         mcmc::mcmc_settings settings;
@@ -255,7 +272,10 @@ void dsge_gensys_R::estim_mcmc_R(const arma::vec& initial_vals, int n_pop, int n
     }
 }
 
-SEXP dsge_gensys_R::IRF_R(int n_irf_periods)
+//
+
+SEXP
+dsge_gensys_R::IRF_R(int n_irf_periods)
 {
     try {
         arma::cube irf_vals = this->IRF(n_irf_periods);
diff --git a/src/modules/dsgevar_gensys_R.cpp b/src/modules/dsgevar_gensys_R.cpp
index c3eaf94..25fab4f 100644
--- a/src/modules/dsgevar_gensys_R.cpp
+++ b/src/modules/dsgevar_gensys_R.cpp
@@ -61,18 +61,15 @@ RCPP_MODULE(dsgevar_gensys_module)
         .field( "mcmc_initial_ub", &dsgevar_gensys_R::mcmc_initial_ub )
 
         .property( "dsge_draws", &dsgevar_gensys_R::get_dsge_draws )
+        .property( "lrem", &dsgevar_gensys_R::get_lrem_R, &dsgevar_gensys_R::set_lrem_R )
+        .property( "dsge", &dsgevar_gensys_R::get_dsge_R, &dsgevar_gensys_R::set_dsge_R )
 
         .method( "build", &dsgevar_gensys_R::build_R )
 
+        .method( "get_model_fn", &dsgevar_gensys_R::get_model_fn )
         .method( "set_model_fn", &dsgevar_gensys_R::set_model_fn )
         .method( "eval_model", &dsgevar_gensys_R::eval_model )
 
-        .method( "get_lrem", &dsgevar_gensys_R::get_lrem_R )
-        .method( "set_lrem", &dsgevar_gensys_R::set_lrem_R )
-
-        .method( "get_dsge", &dsgevar_gensys_R::get_dsge_R )
-        .method( "set_dsge", &dsgevar_gensys_R::set_dsge_R )
-
         .method( "set_bounds", &dsgevar_gensys_R::set_bounds_R )
         .method( "set_prior", &dsgevar_gensys_R::set_prior_R )
 
@@ -86,22 +83,10 @@ RCPP_MODULE(dsgevar_gensys_module)
 //
 // wrapper functions to catch errors and handle memory pointers
 
-void
-dsgevar_gensys_R::build_R(const arma::mat& data_raw, bool cons_term_inp, int p_inp, double lambda_inp)
-{
-    try {
-        this->build(data_raw,cons_term_inp,p_inp,lambda_inp);
-    } catch( std::exception &ex ) {
-        forward_exception_to_r( ex );
-    } catch(...) {
-        ::Rf_error( "BMR: C++ exception (unknown reason)" );
-    }
-}
-
-arma::mat
-dsgevar_gensys_R::get_dsge_draws()
+SEXP
+dsgevar_gensys_R::get_model_fn()
 {
-    return dsge_obj.dsge_draws;
+    return model_fn_SEXP;
 }
 
 void
@@ -109,7 +94,6 @@ dsgevar_gensys_R::set_model_fn(SEXP model_fn_inp)
 {
     try {
         model_fn_SEXP = model_fn_inp;
-        // Rcpp::Function pars_fn = Rcpp::as<Rcpp::Function>(model_fn_SEXP);
         
         dsge_obj.model_fn = std::bind(&dsgevar_gensys_R::model_fn_R,this,std::placeholders::_1,std::placeholders::_2,std::placeholders::_3,std::placeholders::_4,std::placeholders::_5,std::placeholders::_6);
     } catch( std::exception &ex ) {
@@ -119,23 +103,7 @@ dsgevar_gensys_R::set_model_fn(SEXP model_fn_inp)
     }
 }
 
-gensys_R dsgevar_gensys_R::get_lrem_R()
-{
-    gensys_R lrem_obj_out = static_cast<gensys_R&>(dsge_obj.lrem_obj);
-
-    return lrem_obj_out;
-}
-
-void dsgevar_gensys_R::set_lrem_R(gensys_R lrem_obj_inp)
-{
-    try {
-        dsge_obj.lrem_obj = static_cast<bm::gensys&>(lrem_obj_inp);
-    } catch( std::exception &ex ) {
-        forward_exception_to_r( ex );
-    } catch(...) {
-        ::Rf_error( "BMR: C++ exception (unknown reason)" );
-    }
-}
+//
 
 void
 dsgevar_gensys_R::eval_model(Rcpp::NumericVector pars_inp)
@@ -191,14 +159,40 @@ dsgevar_gensys_R::model_fn_R(const arma::vec& pars_inp, bm::gensys& lrem_obj_inp
     }
 }
 
-dsge_gensys_R dsgevar_gensys_R::get_dsge_R()
+// get and set LREM
+
+gensys_R
+dsgevar_gensys_R::get_lrem_R()
+{
+    gensys_R lrem_obj_out = static_cast<gensys_R&>(dsge_obj.lrem_obj);
+
+    return lrem_obj_out;
+}
+
+void
+dsgevar_gensys_R::set_lrem_R(gensys_R lrem_obj_inp)
+{
+    try {
+        dsge_obj.lrem_obj = static_cast<bm::gensys&>(lrem_obj_inp);
+    } catch( std::exception &ex ) {
+        forward_exception_to_r( ex );
+    } catch(...) {
+        ::Rf_error( "BMR: C++ exception (unknown reason)" );
+    }
+}
+
+// get and set DSGE model
+
+dsge_gensys_R
+dsgevar_gensys_R::get_dsge_R()
 {
     dsge_gensys_R dsge_obj_out = static_cast<dsge_gensys_R&>(dsge_obj);
 
     return dsge_obj_out;
 }
 
-void dsgevar_gensys_R::set_dsge_R(dsge_gensys_R dsge_obj_inp)
+void
+dsgevar_gensys_R::set_dsge_R(dsge_gensys_R dsge_obj_inp)
 {
     try {
         dsge_obj = static_cast<bm::dsge<bm::gensys>&>(dsge_obj_inp);
@@ -209,9 +203,30 @@ void dsgevar_gensys_R::set_dsge_R(dsge_gensys_R dsge_obj_inp)
     }
 }
 
+// build and access
+
+void
+dsgevar_gensys_R::build_R(const arma::mat& data_raw, bool cons_term_inp, int p_inp, double lambda_inp)
+{
+    try {
+        this->build(data_raw,cons_term_inp,p_inp,lambda_inp);
+    } catch( std::exception &ex ) {
+        forward_exception_to_r( ex );
+    } catch(...) {
+        ::Rf_error( "BMR: C++ exception (unknown reason)" );
+    }
+}
+
+arma::mat
+dsgevar_gensys_R::get_dsge_draws()
+{
+    return dsge_obj.dsge_draws;
+}
+
 // set bounds and prior
 
-void dsgevar_gensys_R::set_bounds_R(arma::vec lower_bounds_inp, arma::vec upper_bounds_inp)
+void
+dsgevar_gensys_R::set_bounds_R(arma::vec lower_bounds_inp, arma::vec upper_bounds_inp)
 {
     try {
         const int n_vals = lower_bounds_inp.n_elem;
@@ -234,7 +249,8 @@ void dsgevar_gensys_R::set_bounds_R(arma::vec lower_bounds_inp, arma::vec upper_
     }
 }
 
-void dsgevar_gensys_R::set_prior_R(const arma::uvec& prior_form_inp, const arma::mat& prior_pars_inp)
+void
+dsgevar_gensys_R::set_prior_R(const arma::uvec& prior_form_inp, const arma::mat& prior_pars_inp)
 {
     try {
         this->set_prior(prior_form_inp,prior_pars_inp);
@@ -247,7 +263,8 @@ void dsgevar_gensys_R::set_prior_R(const arma::uvec& prior_form_inp, const arma:
 
 //
 
-SEXP dsgevar_gensys_R::estim_mode_R(const arma::vec& initial_vals)
+SEXP
+dsgevar_gensys_R::estim_mode_R(const arma::vec& initial_vals)
 {
     try {
         optim::opt_settings settings;
@@ -271,7 +288,8 @@ SEXP dsgevar_gensys_R::estim_mode_R(const arma::vec& initial_vals)
     return R_NilValue;
 }
 
-void dsgevar_gensys_R::estim_mcmc_R(const arma::vec& initial_vals, int n_pop, int n_gen, int n_burnin)
+void
+dsgevar_gensys_R::estim_mcmc_R(const arma::vec& initial_vals, int n_pop, int n_gen, int n_burnin)
 {
     try {
         mcmc::mcmc_settings settings;
@@ -295,7 +313,10 @@ void dsgevar_gensys_R::estim_mcmc_R(const arma::vec& initial_vals, int n_pop, in
     }
 }
 
-SEXP dsgevar_gensys_R::IRF_R(int n_irf_periods)
+//
+
+SEXP
+dsgevar_gensys_R::IRF_R(int n_irf_periods)
 {
     try {
         arma::cube irf_vals = this->IRF(n_irf_periods);
diff --git a/tests/dsge/gensys/nkm_dsge_gensys.R b/tests/dsge/gensys/nkm_dsge_gensys.R
new file mode 100644
index 0000000..ec55809
--- /dev/null
+++ b/tests/dsge/gensys/nkm_dsge_gensys.R
@@ -0,0 +1,74 @@
+
+#
+
+rm(list=ls())
+
+library(BMR)
+
+source("nkm_model.R")
+
+#
+
+data(BMRVARData)
+dsgedata <- USMacroData[24:211,-c(1,3)]
+dsgedata <- as.matrix(dsgedata)
+for(i in 1:2){
+    dsgedata[,i] <- dsgedata[,i] - mean(dsgedata[,i])
+}
+
+#
+
+obj <- new(dsge_gensys)
+obj$set_model_fn(nkm_model_simple)
+
+x <- c(1)
+
+obj$eval_model(x)
+
+lrem_obj = obj$lrem
+lrem_obj$solve()
+
+lrem_obj$shocks_cov <- matrix(c(1,0,0,0.125),2,2,byrow=TRUE)
+
+sim_data <- lrem_obj$simulate(200,800)$sim_vals
+
+sim_data <- cbind(sim_data[,3],sim_data[,5])
+
+#
+
+prior_pars <- cbind(c(1.0),
+                    c(0.05))
+
+prior_form <- c(1)
+
+obj$set_prior(prior_form,prior_pars)
+
+#
+
+par_bounds <- cbind(c(-Inf),
+                    c( Inf))
+
+opt_bounds <- cbind(c(0.5),
+                    c(3.0))
+
+obj$set_bounds(opt_bounds[,1],opt_bounds[,2])
+
+obj$opt_initial_lb <- opt_bounds[,1]
+obj$opt_initial_ub <- opt_bounds[,2]
+
+#
+
+obj$estim_data = sim_data;
+
+#obj$estim_mode(x)
+
+obj$mcmc_initial_lb <- opt_bounds[,1]
+obj$mcmc_initial_ub <- opt_bounds[,2]
+
+obj$estim_mcmc(x,50,100,100)
+
+varnames <- c("Output Gap","Output","Inflation","Natural Int","Nominal Int","Labour Supply",
+              "Technology","MonetaryPolicy")
+
+plot(obj,parnames="eta",save=FALSE)
+IRF(obj,20,varnames=varnames,save=FALSE)
diff --git a/tests/dsge/gensys/nkm_dsgevar_gensys.R b/tests/dsge/gensys/nkm_dsgevar_gensys.R
new file mode 100644
index 0000000..27dc665
--- /dev/null
+++ b/tests/dsge/gensys/nkm_dsgevar_gensys.R
@@ -0,0 +1,79 @@
+
+#
+
+rm(list=ls())
+
+library(BMR)
+
+source("nkm_model.R")
+
+#
+
+data(BMRVARData)
+dsgedata <- USMacroData[24:211,-c(1,3)]
+dsgedata <- as.matrix(dsgedata)
+for(i in 1:2){
+    dsgedata[,i] <- dsgedata[,i] - mean(dsgedata[,i])
+}
+
+#
+
+obj <- new(dsgevar_gensys)
+obj$set_model_fn(nkm_model_simple)
+
+x <- c(1)
+
+obj$eval_model(x)
+
+#
+
+lrem_obj = obj$lrem
+lrem_obj$solve()
+
+lrem_obj$shocks_cov <- matrix(c(1,0,0,0.125),2,2,byrow=TRUE)
+
+sim_data <- lrem_obj$simulate(200,800)$sim_vals
+
+sim_data <- cbind(sim_data[,3],sim_data[,5])
+
+#
+
+prior_pars <- cbind(c(1.0),
+                    c(0.05))
+
+prior_form <- c(1)
+
+obj$set_prior(prior_form,prior_pars)
+
+#
+
+par_bounds <- cbind(c(-Inf),
+                    c( Inf))
+
+opt_bounds <- cbind(c(0.5),
+                    c(3.0))
+
+obj$set_bounds(opt_bounds[,1],opt_bounds[,2])
+
+obj$opt_initial_lb <- opt_bounds[,1]
+obj$opt_initial_ub <- opt_bounds[,2]
+
+#
+
+cons_term <- TRUE
+p <- 1
+lambda <- 1.0
+
+obj$build(sim_data,cons_term,p,lambda);
+
+#obj$estim_mode(x)
+
+obj$mcmc_initial_lb <- opt_bounds[,1]
+obj$mcmc_initial_ub <- opt_bounds[,2]
+
+obj$estim_mcmc(x,50,100,100)
+
+plot(obj,parnames="eta",save=FALSE)
+IRF(obj,20,varnames=colnames(dsgedata),save=FALSE)
+
+
diff --git a/tests/dsgevar_gensys_simple.R b/tests/dsge/gensys/nkm_model.R
similarity index 69%
rename from tests/dsgevar_gensys_simple.R
rename to tests/dsge/gensys/nkm_model.R
index 4ac3b97..a583e3d 100644
--- a/tests/dsgevar_gensys_simple.R
+++ b/tests/dsge/gensys/nkm_model.R
@@ -1,41 +1,30 @@
 
 #
-
-rm(list=ls())
-library(BMR)
-
-#
-
-data(BMRVARData)
-dsgedata <- USMacroData[24:211,-c(1,3)]
-dsgedata <- as.matrix(dsgedata)
-for(i in 1:2){
-    dsgedata[,i] <- dsgedata[,i] - mean(dsgedata[,i])
-}
-
-#
-
-nkm_model_fn <- function(parameters){
-    alpha    <- 0.33
-    vartheta <- 6
-    beta     <- 0.99
-    theta    <- 0.6667
+# Simple NK Model
+
+nkm_model_fn <- function(parameters)
+{
+    alpha    <- parameters[1]
+    beta     <- parameters[2]
+    vartheta <- parameters[3]
+    theta    <- parameters[4]
+    
+    eta    <- parameters[5]               
+    phi    <- parameters[6]
+    phi_pi <- parameters[7]
+    phi_y  <- parameters[8]
+    rho_a  <- parameters[9]
+    rho_v  <- parameters[10]
     
-    eta    <- parameters[1]               
-    phi    <- 1                  
-    phi_pi <- 1.5             
-    phi_y  <- 0.5/4
-    rho_a  <- 0.90
-    rho_v  <- 0.5
+    sigma_T <- parameters[11]
+    sigma_M <- parameters[12]
     
     BigTheta <- (1-alpha)/(1-alpha+alpha*vartheta)
     kappa    <- (((1-theta)*(1-beta*theta))/(theta))*BigTheta*((1/eta)+((phi+alpha)/(1-alpha)))
     psi      <- (eta*(1+phi))/(1-alpha+eta*(phi + alpha))
     
-    sigma_T <- 1
-    sigma_M <- 0.25
-    
     #
+    
     G0_47 <- (1/eta)*psi*(rho_a - 1)
     #Order:                 yg,       y,      pi,      rn,       i,       n,       a,       v,  yg_t+1,  pi_t+1
     Gamma0 <- rbind(c(      -1,       0,       0,     eta,  -eta/4,       0,       0,       0,       1,   eta/4),
@@ -85,60 +74,28 @@ nkm_model_fn <- function(parameters){
     return=list(Gamma0=Gamma0,Gamma1=Gamma1,GammaC=C,Psi=Psi,Pi=Pi,shocks_cov_out=shocks,C_out=ObsCons,H_out=ObserveMat,R_out=MeasErrs)
 }
 
-obj <- new(dsgevar_gensys)
-obj$set_model_fn(nkm_model_fn)
-
-x <- c(1)
-
-obj$eval_model(x)
-
-lrem_obj = obj$get_lrem()
-lrem_obj$solve()
-
-lrem_obj$shocks_cov <- matrix(c(1,0,0,0.125),2,2,byrow=TRUE)
-
-sim_data <- lrem_obj$simulate(200,800)$sim_vals
-
-sim_data <- cbind(sim_data[,3],sim_data[,5])
-
-#
-
-prior_pars <- cbind(c(1.0),
-                    c(0.05))
-
-prior_form <- c(1)
-
-obj$set_prior(prior_form,prior_pars)
-
-#
-
-par_bounds <- cbind(c(-Inf),
-                    c( Inf))
-
-opt_bounds <- cbind(c(0.5),
-                    c(3.0))
-
-obj$set_bounds(opt_bounds[,1],opt_bounds[,2])
-
-obj$opt_initial_lb <- opt_bounds[,1]
-obj$opt_initial_ub <- opt_bounds[,2]
-
-#
-
-cons_term <- FALSE
-p <- 1
-lambda <- 1.0
-
-obj$build(sim_data,cons_term,p,lambda);
-
-#obj$estim_mode(x)
-
-obj$mcmc_initial_lb <- opt_bounds[,1]
-obj$mcmc_initial_ub <- opt_bounds[,2]
-
-obj$estim_mcmc(x,50,100,100)
-
-plot(obj,parnames="eta",save=FALSE)
-IRF(obj,20,varnames=colnames(dsgedata),save=FALSE)
-
-
+nkm_model_simple <- function(parameters)
+{
+    eta <- parameters[1]
+    
+    #
+    
+    pars_inp <- numeric(12)
+    
+    pars_inp[1]  <- 0.33   # alpha
+    pars_inp[2]  <- 0.99   # beta
+    pars_inp[3]  <- 6.0    # vartheta
+    pars_inp[4]  <- 0.6667 # theta
+    
+    pars_inp[5]  <- eta
+    pars_inp[6]  <- 1      # phi
+    pars_inp[7]  <- 1.5    # phi_pi
+    pars_inp[8]  <- 0.125  # phi_y
+    pars_inp[9]  <- 0.90   # rho_a
+    pars_inp[10] <- 0.5    # rho_v
+    
+    pars_inp[11] <- 1      # sigma_T
+    pars_inp[12] <- 0.25   # sigma_M
+    
+    return(nkm_model_fn(pars_inp))
+}
diff --git a/tests/dsge_gensys_simple.R b/tests/dsge_gensys_simple.R
deleted file mode 100644
index e378d84..0000000
--- a/tests/dsge_gensys_simple.R
+++ /dev/null
@@ -1,141 +0,0 @@
-
-#
-
-rm(list=ls())
-library(BMR)
-
-#
-
-data(BMRVARData)
-dsgedata <- USMacroData[24:211,-c(1,3)]
-dsgedata <- as.matrix(dsgedata)
-for(i in 1:2){
-    dsgedata[,i] <- dsgedata[,i] - mean(dsgedata[,i])
-}
-
-#
-
-nkm_model_fn <- function(parameters){
-    alpha    <- 0.33
-    vartheta <- 6
-    beta     <- 0.99
-    theta    <- 0.6667
-    
-    eta    <- parameters[1]               
-    phi    <- 1                  
-    phi_pi <- 1.5             
-    phi_y  <- 0.5/4
-    rho_a  <- 0.90
-    rho_v  <- 0.5
-    
-    BigTheta <- (1-alpha)/(1-alpha+alpha*vartheta)
-    kappa    <- (((1-theta)*(1-beta*theta))/(theta))*BigTheta*((1/eta)+((phi+alpha)/(1-alpha)))
-    psi      <- (eta*(1+phi))/(1-alpha+eta*(phi + alpha))
-    
-    sigma_T <- 1
-    sigma_M <- 0.25
-    
-    #
-    G0_47 <- (1/eta)*psi*(rho_a - 1)
-    #Order:                 yg,       y,      pi,      rn,       i,       n,       a,       v,  yg_t+1,  pi_t+1
-    Gamma0 <- rbind(c(      -1,       0,       0,     eta,  -eta/4,       0,       0,       0,       1,   eta/4),
-                    c(   kappa,       0,    -1/4,       0,       0,       0,       0,       0,       0,  beta/4),
-                    c(   phi_y,       0,phi_pi/4,       0,    -1/4,       0,       0,       1,       0,       0),
-                    c(       0,       0,       0,      -1,       0,       0,   G0_47,       0,       0,       0),
-                    c(       0,      -1,       0,       0,       0, 1-alpha,       1,       0,       0,       0),
-                    c(      -1,       1,       0,       0,       0,       0,    -psi,       0,       0,       0),
-                    c(       0,       0,       0,       0,       0,       0,       1,       0,       0,       0),
-                    c(       0,       0,       0,       0,       0,       0,       0,       1,       0,       0),
-                    c(       1,       0,       0,       0,       0,       0,       0,       0,       0,       0),
-                    c(       0,       0,       1,       0,       0,       0,       0,       0,       0,       0))
-    #
-    Gamma1 <- rbind(c(       0,       0,       0,       0,       0,       0,       0,       0,       0,       0),
-                    c(       0,       0,       0,       0,       0,       0,       0,       0,       0,       0),
-                    c(       0,       0,       0,       0,       0,       0,       0,       0,       0,       0),
-                    c(       0,       0,       0,       0,       0,       0,       0,       0,       0,       0),
-                    c(       0,       0,       0,       0,       0,       0,       0,       0,       0,       0),
-                    c(       0,       0,       0,       0,       0,       0,       0,       0,       0,       0),
-                    c(       0,       0,       0,       0,       0,       0,   rho_a,       0,       0,       0),
-                    c(       0,       0,       0,       0,       0,       0,       0,   rho_v,       0,       0),
-                    c(       0,       0,       0,       0,       0,       0,       0,       0,       1,       0),
-                    c(       0,       0,       0,       0,       0,       0,       0,       0,       0,       1))
-    #
-    C <- matrix(0,10,1)
-    #
-    Psi <- matrix(0,10,2)
-    Psi[7,1] <- 1
-    Psi[8,2] <- 1
-    #     
-    Pi <- matrix(0,10,2)
-    Pi[9,1] <- 1
-    Pi[10,2] <- 1
-    #
-    
-    shocks <- matrix(c(sigma_T^2,      0,  
-                       0, sigma_M^2),nrow=2)
-    
-    #
-    
-    ObserveMat <- cbind(c(0,0,1,0,0,0,0,0),
-                        c(0,0,0,0,1,0,0,0))
-    
-    ObsCons  <- matrix(0,nrow=2,ncol=1)
-    MeasErrs <- matrix(0,nrow=2,ncol=2)
-    #
-    return=list(Gamma0=Gamma0,Gamma1=Gamma1,GammaC=C,Psi=Psi,Pi=Pi,shocks_cov_out=shocks,C_out=ObsCons,H_out=ObserveMat,R_out=MeasErrs)
-}
-
-obj <- new(dsge_gensys)
-obj$set_model_fn(nkm_model_fn)
-
-x <- c(1)
-
-obj$eval_model(x)
-
-lrem_obj = obj$get_lrem()
-lrem_obj$solve()
-
-lrem_obj$shocks_cov <- matrix(c(1,0,0,0.125),2,2,byrow=TRUE)
-
-sim_data <- lrem_obj$simulate(200,800)$sim_vals
-
-sim_data <- cbind(sim_data[,3],sim_data[,5])
-
-#
-
-prior_pars <- cbind(c(1.0),
-                    c(0.05))
-
-prior_form <- c(1)
-
-obj$set_prior(prior_form,prior_pars)
-
-#
-
-par_bounds <- cbind(c(-Inf),
-                    c( Inf))
-
-opt_bounds <- cbind(c(0.5),
-                    c(3.0))
-
-obj$set_bounds(opt_bounds[,1],opt_bounds[,2])
-
-obj$opt_initial_lb <- opt_bounds[,1]
-obj$opt_initial_ub <- opt_bounds[,2]
-
-#
-
-obj$estim_data = sim_data;
-
-#obj$estim_mode(x)
-
-obj$mcmc_initial_lb <- opt_bounds[,1]
-obj$mcmc_initial_ub <- opt_bounds[,2]
-
-obj$estim_mcmc(x,50,100,100)
-
-varnames <- c("Output Gap","Output","Inflation","Natural Int","Nominal Int","Labour Supply",
-              "Technology","MonetaryPolicy")
-
-plot(obj,parnames="eta",save=FALSE)
-IRF(obj,20,varnames=varnames,save=FALSE)
diff --git a/tests/bvarm.R b/tests/var/bvarm.R
similarity index 100%
rename from tests/bvarm.R
rename to tests/var/bvarm.R
diff --git a/tests/bvars.R b/tests/var/bvars.R
similarity index 100%
rename from tests/bvars.R
rename to tests/var/bvars.R
diff --git a/tests/bvartvp.R b/tests/var/bvartvp.R
similarity index 100%
rename from tests/bvartvp.R
rename to tests/var/bvartvp.R
diff --git a/tests/bvarw.R b/tests/var/bvarw.R
similarity index 100%
rename from tests/bvarw.R
rename to tests/var/bvarw.R
diff --git a/tests/cvar.R b/tests/var/cvar.R
similarity index 100%
rename from tests/cvar.R
rename to tests/var/cvar.R