diff --git a/DESCRIPTION b/DESCRIPTION
index a90139af341ec0a9a35fb6a27cb2b634e35f68b4..9d12fffd5c1d6236e04fe449f4e008820ff27a9d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: onlineforecast
 Type: Package
 Title: Forecast Modelling for Online Applications
-Version: 0.9.3
+Version: 0.9.4
 Description: A framework for fitting adaptive forecasting models. Provides a way to use forecasts as input to models, e.g. weather forecasts for energy related forecasting. The models can be fitted recursively and can easily be setup for updating parameters when new data arrives. See the included vignettes, the website <https://onlineforecasting.org> and the paper "Short-term heat load forecasting for single family houses" <doi:10.1016/j.enbuild.2013.04.022>.
 License: GPL-3
 Encoding: UTF-8
diff --git a/R/bspline.R b/R/bspline.R
index 39103ed6961806577a5714b16f33bb4b5818dc8f..24aa9e76c6a2069fa49167b3846b54b53c452412 100644
--- a/R/bspline.R
+++ b/R/bspline.R
@@ -68,7 +68,7 @@
 #' @export
 bspline <- function(X, Boundary.knots = NA, intercept = FALSE, df = NULL, knots = NULL, degree = 3, bknots = NA, periodic = FALSE) {
     # bknots is just a short for Boundary.knots and replace if Boundary.knots are not given.
-    if(is.na(Boundary.knots)){
+    if(is.na(Boundary.knots[1])){
         Boundary.knots <- bknots
     }
     # If a list, then call on each element
@@ -81,9 +81,12 @@ bspline <- function(X, Boundary.knots = NA, intercept = FALSE, df = NULL, knots
         nams(val) <- nams(X)
         return(val)
     }
-    # X is a data.frame or matrix
+    # X must be a data.frame or matrix
+    if(!(class(X) %in% c("data.frame","matrix"))){ stop("X must be a data.frame or matrix") }
     # First find the horizons, they are used in the end
     nms <- nams(X)
+    # All columns must be like "k12"
+    #if(length(grep("^[k|h][[:digit:]]+$", nms)) != length(nms)){ stop("All column names must indicate a horizon, so start with k and then an integer") }
     # Run for each horizon and calculate the basis splines
     L <- lapply(1:ncol(X), function(i) {
       if (is.na(Boundary.knots[1])) {
diff --git a/R/data.list.R b/R/data.list.R
index 291ac436a361b88ca6b3e00a013f3843fa9f5501..d57a7ac185bc342560de0c7ed64fa9d135ff2c7a 100644
--- a/R/data.list.R
+++ b/R/data.list.R
@@ -330,68 +330,73 @@ check.data.list <- function(object){
     # Which is data.frame or matrix?
     dfOrMat <- sapply(D, function(x){ (class(x) %in% c("matrix","data.frame"))[1] })
     # Vectors check
-    vecchecks <- c("ok","NAs","length","class")
     vecseq <- which(!dfOrMat & names(dfOrMat) != "t")
-    Observations <- data.frame(matrix("", nrow=length(vecseq), ncol=length(vecchecks), dimnames=list(names(vecseq),vecchecks)), stringsAsFactors=FALSE)
-    Observations$ok <- "V"
-    #
-    for(i in 1:length(vecseq)){
+    Observations <- NA
+    if(length(vecseq) > 0){
+        cat("Observation vectors:\n")
+        vecchecks <- c("ok","NAs","length","class")
+        Observations <- data.frame(matrix("", nrow=length(vecseq), ncol=length(vecchecks), dimnames=list(names(vecseq),vecchecks)), stringsAsFactors=FALSE)
+        Observations$ok <- "V"
         #
-        nm <- names(vecseq)[i]
-        # NAs
-        NAs <- round(max(sum(is.na(D[nm])) / length(D[nm])))
-        Observations$NAs[i] <- pst(NAs,"%")
-        # Check the length
-        if(length(D[[nm]]) != length(D$t)){
-            Observations$length[i] <- length(D[[nm]])
-        }
-        # Its class
-        Observations$class[i] <- class(D[[nm]])
-        # Not ok?
-        if(sum(Observations[i, 3] == "") < 1){
-            Observations$ok[i] <- ""
+        for(i in 1:length(vecseq)){
+            #
+            nm <- names(vecseq)[i]
+            # NAs
+            NAs <- round(max(sum(is.na(D[nm])) / length(D[nm])))
+            Observations$NAs[i] <- pst(NAs,"%")
+            # Check the length
+            if(length(D[[nm]]) != length(D$t)){
+                Observations$length[i] <- length(D[[nm]])
+            }
+            # Its class
+            Observations$class[i] <- class(D[[nm]])
+            # Not ok?
+            if(sum(Observations[i, 3] == "") < 1){
+                Observations$ok[i] <- ""
+            }
         }
+        print(Observations)
     }
     #
     # For forecasts
     dfseq <- which(dfOrMat)
-    dfchecks <- c("ok","maxNAs","meanNAs","nrow","colnames","sameclass","class")
-    Forecasts <- data.frame(matrix("", nrow=length(dfseq), ncol=length(dfchecks), dimnames=list(names(dfseq),dfchecks)), stringsAsFactors=FALSE)
-    Forecasts$ok <- "V"
-    #
-    for(i in 1:length(dfseq)){
+    Forecasts <- NA
+    if(length(dfseq) > 0){
+        cat("\nForecast data.frames or matrices:\n")
+        dfchecks <- c("ok","maxNAs","meanNAs","nrow","colnames","sameclass","class")
+        Forecasts <- data.frame(matrix("", nrow=length(dfseq), ncol=length(dfchecks), dimnames=list(names(dfseq),dfchecks)), stringsAsFactors=FALSE)
+        Forecasts$ok <- "V"
         #
-        nm <- names(dfseq)[i]
-        colnms <- nams(D[[nm]])
-        # max NAs
-        maxNAs <- round(max(sapply(colnms, function(colnm){ 100*sum(is.na(D[[nm]][ ,colnm])) / nrow(D[[nm]]) })))
-        Forecasts$maxNAs[i] <- pst(maxNAs,"%")
-        # Mean NAs
-        meanNAs <- round(mean(sapply(colnms, function(colnm){ 100*sum(is.na(D[[nm]][ ,colnm])) / nrow(D[[nm]]) })))
-        Forecasts$meanNAs[i] <- pst(meanNAs,"%")
-        # Check the number of rows
-        if(nrow(D[[nm]]) != length(D$t)){
-            Forecasts$nrow[i] <- nrow(D[[nm]])
-        }
-        # Check the colnames, are they unique and all k+integer?
-        if(!length(unique(grep("^k[[:digit:]]+$",colnms,value=TRUE))) == length(colnms)){
-            Forecasts$colnames[i] <- "X"
-        }
-        if(!length(unique(sapply(colnms, function(colnm){ class(D[[nm]][ ,colnm]) }))) == 1){
-            Forecasts$sameclass[i] <- "X"
-        }else{
-            Forecasts$class[i] <- class(D[[nm]][ ,1])
-        }
-        # Not ok?
-        if(sum(Forecasts[i, ] == "") < (length(dfchecks)-4)){
-            Forecasts$ok[i] <- ""
+        for(i in 1:length(dfseq)){
+            #
+            nm <- names(dfseq)[i]
+            colnms <- nams(D[[nm]])
+            # max NAs
+            maxNAs <- round(max(sapply(colnms, function(colnm){ 100*sum(is.na(D[[nm]][ ,colnm])) / nrow(D[[nm]]) })))
+            Forecasts$maxNAs[i] <- pst(maxNAs,"%")
+            # Mean NAs
+            meanNAs <- round(mean(sapply(colnms, function(colnm){ 100*sum(is.na(D[[nm]][ ,colnm])) / nrow(D[[nm]]) })))
+            Forecasts$meanNAs[i] <- pst(meanNAs,"%")
+            # Check the number of rows
+            if(nrow(D[[nm]]) != length(D$t)){
+                Forecasts$nrow[i] <- nrow(D[[nm]])
+            }
+            # Check the colnames, are they unique and all k+integer?
+            if(!length(unique(grep("^k[[:digit:]]+$",colnms,value=TRUE))) == length(colnms)){
+                Forecasts$colnames[i] <- "X"
+            }
+            if(!length(unique(sapply(colnms, function(colnm){ class(D[[nm]][ ,colnm]) }))) == 1){
+                Forecasts$sameclass[i] <- "X"
+            }else{
+                Forecasts$class[i] <- class(D[[nm]][ ,1])
+            }
+            # Not ok?
+            if(sum(Forecasts[i, ] == "") < (length(dfchecks)-4)){
+                Forecasts$ok[i] <- ""
+            }
         }
+        print(Forecasts)
     }
-    #
-    message("Observation vectors:")
-    print(Observations)
-    message("\nForecast data.frames or matrices:")
-    print(Forecasts)
 
     invisible(list(Observations=Observations, Forecasts=Forecasts))
 }
diff --git a/R/lagdl.R b/R/lagdl.R
new file mode 100644
index 0000000000000000000000000000000000000000..5f6696fe6bca8e303b8c0e8ecaea012afeef2579
--- /dev/null
+++ b/R/lagdl.R
@@ -0,0 +1,28 @@
+## Do this in a separate file to see the generated help:
+#library(devtools)
+#document()
+#load_all(as.package("../../onlineforecast"))
+#?lagdl
+
+#' Lagging by shifting the values back or fourth always returning a data.list.
+#'
+#' This function lags (shifts) the values of the vector. A data.list is always returned with each data.frame lagged with \code{lagdf}.
+#'
+#' 
+#' @title Lagging which returns a data.list
+#' @param x The data.list to be lagged.
+#' @param lagseq The integer(s) setting the lag steps.
+#' @return A data.list.
+#' @seealso \code{\link{lagdl.data.frame}} which is run when \code{x} is a data.frame.
+#' @examples
+#' # The values are simply shifted in each data.frame with lagdf
+#'
+#'@export
+
+lagdl <- function(D, lagseq){
+    iseq <- which(sapply(D,class) %in% c("data.frame","matrix"))
+    D[iseq] <- lapply(iseq, function(i){
+        lagdf(D[[i]], lagseq)
+    })
+    return(D)
+}
diff --git a/R/make_input.R b/R/make_input.R
index ff525d114e1b2d1c60182314a3340b39c9cb1c6c..b7ea05c30ee2a7f794ae0dd9548d8dd46476a0d5 100644
--- a/R/make_input.R
+++ b/R/make_input.R
@@ -28,7 +28,7 @@ make_input <- function(observations, kseq){
     val <- sapply(kseq, function(k){
         observations
     })
-    ## set row and column names
+    # set row and column names
     nams(val) <- paste0('k', kseq)
     return( as.data.frame(val) )
 }
diff --git a/inst/CITATION b/inst/CITATION
index a7bf2f42f7926cf58a1e068d906609aef6040430..38f09800d55540f4c9a1a71b9831d95720626c87 100644
--- a/inst/CITATION
+++ b/inst/CITATION
@@ -15,3 +15,12 @@ citEntry(
       "We are in process of writing a journal paper about the package, but for now we referer to the paper 'Short-term heat load forecasting for single family houses', in which the implemented modelling is described."
   )
 )
+
+citEntry(
+  entry    = "Manual",
+  title    = "{{onlineforecast}: Forecast Modelling for Online Applications}",
+  author   = "Peder Bacher and Hjörleifur G. Bergsteinsson",
+  year     = "2020",
+  note     = "R package version 0.9.3",
+  url      = "https://onlineforecasting.org"
+)
diff --git a/make.R b/make.R
index 25614e49c15013cb506efa20af15ca5e36dff562..cf5912c5c76f446c79c0045b03ed639307a34ede 100644
--- a/make.R
+++ b/make.R
@@ -60,10 +60,10 @@ library(roxygen2)
 # ----------------------------------------------------------------
 # Build the package
 document()
-build(".", vignettes=TRUE)
+build(".", vignettes=FALSE)
 
 # Install it
-install.packages("../onlineforecast_0.9.3.tar.gz")
+install.packages("../onlineforecast_0.9.4.tar.gz")
 
 library(onlineforecast)
 # ----------------------------------------------------------------
diff --git a/vignettes/setup-and-use-model.Rmd b/vignettes/setup-and-use-model.Rmd
index 73b537ec5f68ba62b60c6089b5b36bcf3be400fa..2b04bc6329ac3e63af5937f96ad5cc61ba4f829c 100644
--- a/vignettes/setup-and-use-model.Rmd
+++ b/vignettes/setup-and-use-model.Rmd
@@ -362,6 +362,7 @@ There are quite a few functions available for input transformations:
 - `one()` generates an matrix of ones (for including an intercept).
 - `fs()` generate Fourier series for modelling harmonic functions.
 - `bspline()` wraps the `bs()` function for generating base splines.
+- `pbspline()` wraps the `pbs()` function for generating periodic base splines.
 - `AR()` generates auto-regressive model inputs.
 
 and they can even be combined, see more details in [onlineforecasting] and in their help