diff --git a/Project.toml b/Project.toml
index 87872b8..b9c0672 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "StateSpaceEcon"
 uuid = "e4c825b0-b65c-11ea-0b5a-6176b64e7b7f"
 authors = ["Atai Akunov <aakunov@bankofcanada.ca>", "Boyan Bejanov <bbejanov@bankofcanada.ca>", "Nicholas Labelle St-Pierre <labs@bankofcanada.ca>"]
-version = "0.1.0"
+version = "0.1.1"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
@@ -11,17 +11,19 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelBaseEcon = "1d16192e-b65e-11ea-11ed-0789cee22d2f"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 TimeSeriesEcon = "8b6756d2-c55c-11ea-2998-5f67ea17da60"
 
 [compat]
 DiffResults = "1.0"
 ForwardDiff = "0.10"
 ModelBaseEcon = "0.1"
-TimeSeriesEcon = "0.1"
+TimeSeriesEcon = "0.2"
 julia = "1.0"
 
 [extras]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
diff --git a/src/Plans.jl b/src/Plans.jl
index 4a130ef..6d7ed74 100644
--- a/src/Plans.jl
+++ b/src/Plans.jl
@@ -62,6 +62,7 @@ range of the plan is augmented to include periods before and after the given
 range, over which initial and final conditions will be applied.
 
 """
+Plan(model::Model, r::MIT) = Plan(model, r:r)
 function Plan(model::Model, range::AbstractUnitRange)
     if !(eltype(range) <: MIT)
         range = UnitRange{MIT{Unit}}(range)
@@ -82,8 +83,8 @@ Base.axes(p::Plan) = (p.range,)
 Base.length(p::Plan) = length(p.range)
 Base.IndexStyle(::Plan) = IndexLinear()
 
-@inline _offset(p::Plan{T}, idx::T) where T = 1 - first(p.range) + idx
-@inline _offset(p::Plan{T}, idx::AbstractUnitRange{T}) where T = 1 - first(p.range) .+ idx
+@inline _offset(p::Plan{T}, idx::T) where T = idx - first(p.range) + 1
+@inline _offset(p::Plan{T}, idx::AbstractUnitRange{T}) where T = idx .- first(p.range) .+ 1
 
 ## Integer index is taken as is
 Base.getindex(p::Plan, idx::Int) = p[p.range[idx]]
@@ -91,7 +92,7 @@ Base.getindex(p::Plan, idx::AbstractUnitRange{Int}) = p[p.range[idx]]
 Base.getindex(p::Plan, idx::Int, ::Val{:inds}) = findall(p.exogenous[idx,:])
 
 # Index of the frequency type returns the list of exogenous symbols
-Base.getindex(p::Plan{T}, idx::T) where T <: MIT = [keys(p.varsshks)[p[_offset(p,idx),Val(:inds)]]...,]
+Base.getindex(p::Plan{T}, idx::T) where T <: MIT = [keys(p.varsshks)[p[_offset(p, idx),Val(:inds)]]...,]
 function Base.getindex(p::Plan{T}, idx::T, ::Val{:inds}) where T <: MIT
     first(p.range) ≤ idx ≤ last(p.range) || throw(BoundsError(p, idx))
     return p[_offset(p, idx), Val(true)]
@@ -102,7 +103,7 @@ Base.getindex(p::Plan{MIT{Unit}}, rng::AbstractUnitRange{Int}) = p[UnitRange{MIT
 function Base.getindex(p::Plan{T}, rng::AbstractUnitRange{T}) where T <: MIT
     rng.start < p.range.start && throw(BoundsError(p, rng.start))
     rng.stop > p.range.stop && throw(BoundsError(p, rng.stop))
-    return Plan{T}(rng, p.varsshks, p.exogenous[_offset(rng), :])
+    return Plan{T}(rng, p.varsshks, p.exogenous[_offset(p, rng), :])
 end
 
 # A range with a model returns a plan trimmed over that range and extended for initial and final conditions.
@@ -119,25 +120,53 @@ Base.setindex!(p::Plan, x, i...) = error("Cannot assign directly. Use `exogenize
 
 Base.summary(io::IO, p::Plan) = print(io, typeof(p), " with range ", p.range)
 
-function Base.show(io::IO, ::MIME"text/plain", p::Plan)
+#  Temporary fix to override bugs in TimeSeriesEcon
+Base.axes(r::AbstractUnitRange{<:MIT}) = (Base.OneTo(length(r)),)
+Base.axes1(r::AbstractUnitRange{<:MIT}) = Base.OneTo(length(r))
+
+# Used in the show() implementation below
+function collapsed_range(p::Plan{T}) where T <: MIT
+    ret = Pair{Union{T,UnitRange{T}},Vector{Symbol}}[]
+    i1 = first(p.range)
+    i2 = i1
+    val = p[i1]
+    make_key() = i1 == i2 ? i1 : i1:i2
+    for i in p.range[2:end]
+        val_i = p[i]
+        if val_i == val
+            i2 = i
+        else
+            push!(ret, make_key() => val)
+            i1 = i2 = i
+            val = val_i
+        end
+    end
+    push!(ret, make_key() => val)
+end
+
+Base.show(io::IO, ::MIME"text/plain", p::Plan) = show(io, p)
+function Base.show(io::IO, p::Plan)
     # 0) show summary before setting :compact
     summary(io, p)
     isempty(p) && return
     print(io, ":")
     nrow, ncol = displaysize(io)
-    if length(p) <= nrow - 5
-        for r in p.range
-            print(io, "\n  ", r, " => ", p[r])
+    cp = collapsed_range(p)
+    # find the longest string left of "=>" for padding
+    maxl = maximum(length("$k") for (k,v) in cp)
+    if length(cp) <= nrow - 5
+        for (r, v) in cp
+            print(io, "\n  ", lpad("$r", maxl, " "), " => ", v)
         end
     else
         top = div(nrow - 5, 2)
-        bot = length(p.range) - nrow + 6 + top
-        for r in p.range[1:top]
-            print(io, "\n  ", r, " => ", p[r])
+        bot = length(cp) - nrow + 6 + top
+        for (r, v) in cp[1:top]
+            print(io, "\n  ", lpad("$r", maxl, " "), " => ", v)
         end
         print(io, "\n   ⋮")
-        for r in p.range[bot:end]
-            print(io, "\n  ", r, " => ", p[r])
+        for (r, v) in cp[bot:end]
+            print(io, "\n  ", lpad("$r", maxl, " "), " => ", v)
         end
     end
 end
@@ -156,7 +185,7 @@ Modify the status of the given variable(s) on the given date(s). If `value` is
 """
 function setplanvalue!(p::Plan{T}, val::Bool, vars::Array{Symbol,1}, date::AbstractUnitRange{T}) where T <: MIT
     firstindex(p) ≤ first(date) && last(date) ≤ lastindex(p) || throw(BoundsError(p, date))
-    idx1 = 1 - firstindex(p) .+ date
+    idx1 = date .- firstindex(p) .+ 1
     for v in vars
         idx2 = get(p.varsshks, v, 0)
         if idx2 == 0
@@ -168,8 +197,9 @@ function setplanvalue!(p::Plan{T}, val::Bool, vars::Array{Symbol,1}, date::Abstr
 end
 setplanvalue!(p::Plan{MIT{Unit}}, val::Bool, vars::Array{Symbol,1}, date::AbstractUnitRange{Int}) = setplanvalue!(p, val, vars, UnitRange{MIT{Unit}}(date))
 setplanvalue!(p::Plan{MIT{Unit}}, val::Bool, vars::Array{Symbol,1}, date::Integer) = setplanvalue!(p, val, vars, ii(date):ii(date))
+setplanvalue!(p::Plan{MIT{Unit}}, val::Bool, vars::Array{Symbol,1}, date::MIT{Unit}) = setplanvalue!(p, val, vars, date:date)
 setplanvalue!(p::Plan{T}, val::Bool, vars::Array{Symbol,1}, date::T) where T <: MIT = setplanvalue!(p, val, vars, date:date)
-setplanvalue!(p::Plan, val::Bool, vars::Array{Symbol,1}, date) = (foreach(d->setplanvalue!(p, val, vars, d), date); p)
+setplanvalue!(p::Plan, val::Bool, vars::Array{Symbol,1}, date) = (foreach(d -> setplanvalue!(p, val, vars, d), date); p)
 
 
 """
@@ -248,6 +278,11 @@ end
 #######################################
 # The user interface to prepare data for simulation.
 
+TimeSeriesEcon.frequencyof(p::Plan) = frequencyof(p.range)
+TimeSeriesEcon.firstdate(p::Plan) = first(p.range)
+TimeSeriesEcon.lastdate(p::Plan) = last(p.range)
+import ..SimData
+
 """
     zeroarray(model, plan)
     zeroarray(model, range)
@@ -296,7 +331,7 @@ See also: [`zeroarray`](@ref), [`zerodict`](@ref), [`steadystatearray`](@ref),
 [`steadystatedict`](@ref)
 
 """
-zerodata(::Model, p::Plan) = NamedTuple{keys(p.varsshks)}(((TSeries(p.range, 0.0) for _ in p.varsshks)...,))
+zerodata(::Model, p::Plan) = SimData(firstdate(p), keys(p.varsshks), zeros(length(p), length(p.varsshks)))
 zerodata(m::Model, rng::AbstractUnitRange) = zerodata(m, Plan(m, rng))
 
 """
@@ -337,7 +372,7 @@ steadystatedict(m::Model, p::Plan) = Dict(string(v) => TSeries(p.range, i <= Mod
     steadystatedata(model, plan)
     steadystatedata(model, range)
 
-Create a dictionary containing a [`TSeries`](@ref) of the appropriate range for each
+Create a [`SimData`](@ref) containing a [`TSeries`](@ref) of the appropriate range for each
 variable in the model for a simulation with the given plan or over the given
 range. The matrix is initialized with the steady state level of each variable.
 If a range is given rather than a plan, it is augmented with periods before and
@@ -348,7 +383,10 @@ See also: [`zeroarray`](@ref), [`zerodict`](@ref), [`steadystatearray`](@ref),
 
 """
 steadystatedata(m::Model, rng::AbstractUnitRange) = steadystatedict(m, Plan(m, rng))
-steadystatedata(m::Model, p::Plan) = NamedTuple{keys(p.varsshks)}(((TSeries(p.range, i ≤ ModelBaseEcon.nvariables(m) ? m.sstate[v] : 0) for (v, i) in pairs(p.varsshks))...,))
+steadystatedata(m::Model, p::Plan) = hcat(
+    SimData(firstdate(p), (), zeros(length(p), 0)); 
+    (v => (i <= ModelBaseEcon.nvariables(m) ? m.sstate[v] : 0) for (v, i) in pairs(p.varsshks))...
+)
 
 #######################################
 # The internal interface to simulations code.
diff --git a/src/StateSpaceEcon.jl b/src/StateSpaceEcon.jl
index e5691d9..e00ad00 100644
--- a/src/StateSpaceEcon.jl
+++ b/src/StateSpaceEcon.jl
@@ -11,6 +11,7 @@ using ModelBaseEcon
 # using ModelBaseEcon.OptionsMod
 # using ModelBaseEcon.Timer
 
+include("simdata.jl")
 include("SteadyStateSolver.jl")
 include("Plans.jl")
 include("misc.jl")
diff --git a/src/simdata.jl b/src/simdata.jl
new file mode 100644
index 0000000..53cacab
--- /dev/null
+++ b/src/simdata.jl
@@ -0,0 +1,354 @@
+
+export SimData
+
+"""
+    SimData
+
+Data structure containing the time series data for a simulation.
+
+It is a collection of [`TSeries`](@ref) of the same frequency and containing
+data for the same range. When used for simulation, the range must include the
+initial conditions, the simulation range and the final conditions, although it
+could extend beyond that. It must contain time series for all variables and
+shocks in the model, although it might contain other time series.
+
+"""
+struct SimData{F <: Frequency,N <: NamedTuple,C <: AbstractMatrix{Float64}} <: AbstractMatrix{Float64}
+    firstdate::MIT{F}
+    columns::N
+    values::C
+
+    # inner constructor enforces constraints.
+    function SimData(fd, names, values)
+        if length(names) != size(values, 2)
+            throw(ArgumentError("Number of names and columns don't match: $(length(names)) ≠ $(size(values, 2))."))
+        end
+        columns = NamedTuple{names}([TSeries(fd, view(values, :, i)) for i in 1:length(names) ])
+        new{frequencyof(fd),typeof(columns),typeof(values)}(fd, columns, values)
+    end
+end
+
+# External constructors
+# Coming soon stay tuned
+
+TimeSeriesEcon.firstdate(sd::SimData) = getfield(sd, :firstdate)
+TimeSeriesEcon.lastdate(sd::SimData) = (firstdate(sd) - 1) + size(_values(sd), 1) 
+TimeSeriesEcon.frequencyof(::SimData{F}) where F <: Frequency = F
+TimeSeriesEcon.mitrange(sd::SimData) = (firstdate(sd) - 1) .+ Base.axes1(sd)
+
+@inline _columns(sd::SimData) = getfield(sd, :columns)
+@inline _names(sd::SimData) = keys(getfield(sd, :columns))
+@inline _values(sd::SimData) = getfield(sd, :values)
+@inline _values_view(sd::SimData, i1, i2) = view(getfield(sd, :values), i1, i2)
+@inline _values_slice(sd::SimData, i1, i2) = getindex(getfield(sd, :values), i1, i2)
+
+Base.size(sd::SimData) = size(_values(sd))
+Base.IndexStyle(sd::SimData) = IndexStyle(_values(sd))
+Base.getindex(sd::SimData, i1...) = getindex(_values(sd), i1...)
+Base.setindex!(sd::SimData, val, i1...) = setindex!(_values(sd), val, i1...)
+
+Base.similar(sd::SimData) = SimData(firstdate(sd), _names(sd), similar(_values(sd)))
+
+Base.:(==)(a::SimData, b::SimData) = frequencyof(a) == frequencyof(b) && firstdate(a) == firstdate(b) && _values(a) == _values(b)
+
+Base.dataids(sd::SimData) = Base.dataids(_values(sd))
+
+# Define dot access to columns
+Base.propertynames(sd::SimData) = _names(sd)
+function Base.getproperty(sd::SimData, col::Symbol)
+    if col in _names(sd)
+        return getfield(_columns(sd), col)
+    else
+        throw(BoundsError(sd, [col,]))
+    end
+end
+
+# Access to columns by [:xyz] notation
+Base.getindex(sd::SimData, col::Symbol) = getproperty(sd, col)
+Base.getindex(sd::SimData, col::AbstractString) = getproperty(sd, Symbol(col))
+
+function Base.setproperty!(sd::SimData,  col::Symbol, val)
+    try
+        col = getproperty(sd, col)
+    catch BoundsError
+        error("Cannot assign new data column this way. Use hcat(sd, column=value, ...) instead")
+    end
+    if Base.mightalias(col, val)
+        val = copy(val)
+    end
+    setindex!(col, val, :)
+end
+
+Base.getindex(sd::SimData, cols::AbstractArray{Symbol}) = getindex(sd, tuple(cols...))
+function Base.getindex(sd::SimData, cols::NTuple{N,Symbol}) where N
+    SimData(firstdate(sd), tuple(cols...), hcat((getproperty(sd, c) for c in cols)...))
+end
+
+function Base.hcat(sd::SimData; KW...)
+    l1 = size(sd, 1)
+    as_vect(v::Number) = fill(Float64(v), l1)
+    as_vect(v) = v
+    names = (_names(sd)..., keys(KW)...)
+    vals = hcat(_values(sd), (as_vect(v) for v in values(KW))...)
+    return SimData(firstdate(sd), names, vals)
+end
+
+# Indexing access with MIT
+
+function check_frequency(a, b)
+    Fa = frequencyof(a)
+    Fb = frequencyof(b)
+    Fa == Fb || throw(ArgumentError("wrong frequency: expected $Fa got $Fb."))
+end
+
+macro if_same_frequency(a, b, expr)
+    :( frequencyof($a) == frequencyof($b) ? $(expr) : throw(ArgumentError("Wrong frequency")) ) |> esc
+end
+
+# A single MIT materializes as a NamedTuple (row of the matrix with column names attached to the values)
+function Base.getindex(sd::SimData, i1::MIT) 
+    check_frequency(sd, i1)
+    if firstdate(sd) <= i1 <= lastdate(sd)
+        return @inbounds NamedTuple{_names(sd)}(_values_view(sd, i1 - firstdate(sd) + 1, :))
+    else
+        throw(BoundsError(sd, i1))
+    end
+end
+
+# Modifying a row in a table -> one must pass in a vector
+function Base.setindex!(sd::SimData, val::AbstractVector{<:Real}, i1::MIT) 
+    check_frequency(sd, i1)
+    if firstdate(sd) <= i1 <= lastdate(sd)
+        row = i1 - firstdate(sd) + 1
+        setindex!(_values_view(sd, row, :), val, :)
+        return _values(sd)[row,:]
+    else
+        throw(BoundsError(sd, i1))
+    end
+end
+
+function Base.setindex!(sd::SimData, val::NamedTuple, i1::MIT)
+    check_frequency(sd, i1)
+    if firstdate(sd) <= i1 <= lastdate(sd)
+        for (n, v) in pairs(val)
+            setindex!(getproperty(sd, n), v, i1)
+        end
+        return sd[i1]
+    else
+        throw(BoundsError(sd, i1))
+    end
+    return sd[i1]
+end
+
+# A selection of several rows returns a slice from the original SimData
+function Base.getindex(sd::SimData, i1::AbstractUnitRange{<:MIT})
+    check_frequency(sd, i1)
+    if firstdate(sd) <= minimum(i1) <= maximum(i1) <= lastdate(sd)
+        return SimData(first(i1), _names(sd), _values_slice(sd, i1 .- firstdate(sd) .+ 1, :))
+    else
+        throw(BoundsError(sd, i1))
+    end
+end
+
+function Base.setindex!(sd::SimData, val, i1::AbstractUnitRange{<:MIT})
+    check_frequency(sd, i1)
+    if firstdate(sd) <= minimum(i1) <= maximum(i1) <= lastdate(sd)
+        rows = i1 .- firstdate(sd) .+ 1
+        setindex!(_values_view(sd, rows, :), val, :, :)
+        return sd[rows,:]
+    else
+        throw(BoundsError(sd, i1))
+    end
+end
+
+
+Base.getindex(sd::SimData, r, c::Symbol) = getproperty(sd, c)[r]
+Base.getindex(sd::SimData, r::MIT, c::NTuple{N,Symbol}) where N = (a = sd[r]; NamedTuple{c}([a[cc] for cc in c]))
+Base.getindex(sd::SimData, r::MIT, c::AbstractVector{Symbol}) = (a = sd[r]; NamedTuple{tuple(c...)}([a[cc] for cc in c]))
+Base.getindex(sd::SimData, r::AbstractUnitRange{<:MIT}, c::AbstractVector{Symbol}) = getindex(sd, r, tuple(c...))
+function Base.getindex(sd::SimData, r::AbstractUnitRange{<:MIT}, c::NTuple{N,Symbol}) where N 
+    check_frequency(sd, r)
+    if firstdate(sd) <= first(r) <= last(r) <= lastdate(sd)
+        col_inds = indexin(c, [_names(sd)...])
+        for (cc, cn) in zip(c, col_inds)
+            if cn === nothing
+                throw(BoundsError(sd, [cc]))
+            end
+        end
+        return SimData(first(r), c, _values_slice(sd, r .- firstdate(sd) .+ 1, col_inds))
+    else
+        throw(BoundsError(sd, r))
+    end
+end
+
+Base.setindex!(sd::SimData, v, r, c::Symbol) = setindex!(getproperty(sd, c), v, r)
+function Base.setindex!(sd::SimData, vals, r::MIT, c::NTuple{N,Symbol}) where N 
+    for (cc, vv) in zip(c, vals)
+        setindex!(sd, vv, r, cc)
+    end
+    return sd[r, c]
+end
+Base.setindex!(sd::SimData, vals, r::MIT, c::AbstractVector{Symbol}) = setindex!(sd, vals, r, tuple(c...))
+Base.setindex!(sd::SimData, vals, r::AbstractUnitRange{<:MIT}, c::AbstractVector{Symbol}) = setindex!(sd, vals, r, tuple(c...))
+function Base.setindex!(sd::SimData, vals, r::AbstractUnitRange{<:MIT}, c::NTuple{N,Symbol}) where N
+    check_frequency(sd, r)
+    if firstdate(sd) <= first(r) <= last(r) <= lastdate(sd)
+        cols = indexin(c, [_names(sd)...])
+        for (cc, cn) in zip(c, cols)
+            if cn === nothing
+                throw(BoundsError(sd, [cc]))
+            end
+        end
+        rows = r .- firstdate(sd) .+ 1
+        if vals isa Number
+            vals = fill(Float64(vals), length(r) * length(c))
+        end
+        setindex!(_values_view(sd, rows, cols), vals, :, :)
+    else
+        throw(BoundsError(sd, r))
+    end
+end
+
+#### Pretty printing
+
+sprint_names(names) = length(names) > 10 ? "$(length(names)) variables" : "variables (" * join(names, ",") * ")"
+Base.summary(io::IO, sd::SimData) = isempty(sd) ?
+        print(IOContext(io, :limit => true), "Empty SimData with ", sprint_names(_names(sd)), " in ", mitrange(sd)) :
+        print(IOContext(io, :limit => true), size(sd, 1), '×', size(sd, 2), " SimData with ", sprint_names(_names(sd)), " in ", mitrange(sd))
+
+
+Base.show(io::IO, ::MIME"text/plain", sd::SimData) = show(io, sd)
+function Base.show(io::IO, sd::SimData)
+    summary(io, sd)
+    isempty(sd) && return
+    print(io, ":")
+    limit = get(io, :limit, true)
+    nval, nsym = size(sd)
+
+    from = firstdate(sd)
+    dheight, dwidth = displaysize(io)
+    if get(io, :compact, nothing) === nothing
+        io = IOContext(io, :compact => true)
+    end
+    dwidth -= 11
+
+    names_str = let 
+        names = map(_names(sd)) do n
+            sn = string(n)
+            return sn[1:min(10, length(sn))]
+        end
+        reshape([names...], 1, :)
+    end
+    sd_with_names = [names_str; _values(sd)]
+    A = Base.alignment(io, sd_with_names, axes(sd_with_names, 1), 1:nsym, dwidth, dwidth, 2)
+
+    all_cols = true
+    if length(A) ≠ nsym
+        dwidth = div(dwidth - 1, 2)
+        AL = Base.alignment(io, sd_with_names, axes(sd_with_names, 1), 1:nsym, dwidth, dwidth, 2)
+        AR = reverse(Base.alignment(io, sd_with_names, axes(sd_with_names, 1), reverse(1:nsym), dwidth, dwidth, 2))
+        Linds = [1:length(AL)...]
+        Rinds = [nsym - length(AR) + 1:nsym...]
+        all_cols = false
+    end
+
+    local vdots = "\u22ee"
+    local hdots = " \u2026 "
+    local ddots = " \u22f1 "
+
+    print_aligned_val(io, v, (al, ar), showsep=true; sep=showsep ? "  " : "") = begin
+        sv = sprint(print, v, context=io, sizehint=0)
+        if v isa Number
+            vl, vr = Base.alignment(io, v)
+        else
+            if length(sv) > al + ar
+                sv = sv[1:al + ar - 1] * '…'
+            end
+            vl, vr = al, length(sv) - al
+        end
+        print(io, repeat(" ", al - vl), sv, repeat(" ", ar - vr), sep)
+    end
+
+    print_colnames(io, Lcols, LAligns, Rcols=[], RAligns=[]) = begin
+        local nLcols = length(Lcols)
+        local nRcols = length(Rcols)
+        for (i, (col, align)) in enumerate(zip(Lcols, LAligns))
+            print_aligned_val(io, names[col], align, i < nLcols)
+        end
+        nRcols == 0 && return
+        print(io, hdots)
+        for (i, (col, align)) in enumerate(zip(Rcols, RAligns))
+            print_aligned_val(io, names[col], align, i < nRcols)
+        end
+    end
+
+    print_rows(io, rows, Lcols, LAligns, Rcols=[], RAligns=[]) = begin
+        local nLcols = length(Lcols)
+        local nRcols = length(Rcols)
+        for row in rows
+            mit = from + (row - 1)
+            print(io, '\n', lpad(mit, 8), " : ")
+            for (i, (val, align)) in enumerate(zip(vals[row, Lcols], LAligns))
+                print_aligned_val(io, val, align, i < nLcols)
+            end
+            nRcols == 0 && continue
+            print(io, hdots)
+            for (i, (val, align)) in enumerate(zip(vals[row, Rcols], RAligns))
+                print_aligned_val(io, val, align, i < nRcols)
+            end
+        end
+    end
+
+    print_vdots(io, Lcols, LAligns, Rcols=[], RAligns=[]) = begin
+        print(io, '\n', repeat(" ", 11))
+        local nLcols = length(Lcols)
+        local nRcols = length(Rcols)
+        for (i, (col, align)) in enumerate(zip(Lcols, LAligns))
+            print_aligned_val(io, vdots, align, i < nLcols)
+        end
+        nRcols == 0 && return
+        print(io, ddots)
+        for (i, (col, align)) in enumerate(zip(Rcols, RAligns))
+            print_aligned_val(io, vdots, align, i < nRcols)
+        end
+    end
+
+    names = _names(sd)
+    vals = _values(sd)
+    if !limit
+        print(io, "\n", repeat(" ", 11))
+        print_colnames(io, 1:nsym, A)
+        print_rows(io, 1:nval, 1:nsym, A)
+    elseif nval > dheight - 6 # all rows don't fit
+                # unable to show all rows
+        if all_cols
+            print(io, "\n", repeat(" ", 11))
+            print_colnames(io, 1:nsym, A)
+            top = div(dheight - 6, 2)
+            print_rows(io, 1:top, 1:nsym, A)
+            print_vdots(io, 1:nsym, A)
+            bot = nval - dheight + 7 + top
+            print_rows(io, bot:nval, 1:nsym, A)
+        else # not all_cols
+            print(io, "\n", repeat(" ", 11))
+            print_colnames(io, Linds, AL, Rinds, AR)
+            top = div(dheight - 6, 2)
+            print_rows(io, 1:top, Linds, AL, Rinds, AR)
+            print_vdots(io, Linds, AL, Rinds, AR)
+            bot = nval - dheight + 7 + top
+            print_rows(io, bot:nval, Linds, AL, Rinds, AR)
+        end # all_cols
+    else # all rows fit
+        if all_cols
+            print(io, '\n', repeat(" ", 11))
+            print_colnames(io, 1:nsym, A)
+            print_rows(io, 1:nval, 1:nsym, A)
+        else
+            print(io, '\n', repeat(" ", 11))
+            print_colnames(io, Linds, AL, Rinds, AR)
+            print_rows(io, 1:nval, Linds, AL, Rinds, AR)
+        end
+    end
+end
+
+
diff --git a/src/stackedtime/simulate.jl b/src/stackedtime/simulate.jl
index 3f4ac17..4916b23 100644
--- a/src/stackedtime/simulate.jl
+++ b/src/stackedtime/simulate.jl
@@ -16,7 +16,7 @@ Solve the simulation problem.
 """
 function sim_nr!(x::AbstractArray{Float64}, sd::AbstractSolverData,
                 maxiter::Int64, tol::Float64, verbose::Bool,
-                sparse_solver::Function = (A, b)->A \ b)
+                sparse_solver::Function=(A, b) -> A \ b)
     for it = 1:maxiter
         @timer Fx, Jx = global_RJ(x, x, sd)
         @timer nFx = norm(Fx, Inf)
@@ -86,16 +86,16 @@ Run a simulation for the given model, simulation plan and exogenous data.
 
 """
 function simulate(m::Model, p::Plan, exog_data::AbstractArray{Float64,2};
-                    initial_guess::AbstractArray{Float64,2} = zeros(0, 0),
-                    linearize::Bool = false,
-                    deviation::Bool = false,
-                    anticipate::Bool = true,
-                    verbose::Bool = m.options.verbose,
-                    tol::Float64 = m.options.tol,
-                    maxiter::Int64 = m.options.maxiter,
-                    fctype::FCType = fcgiven,
-                    expectation_horizon::Union{Nothing,Int64} = nothing,
-                    sparse_solver::Function = (A, b)->A \ b
+                    initial_guess::AbstractArray{Float64,2}=zeros(0, 0),
+                    linearize::Bool=false,
+                    deviation::Bool=false,
+                    anticipate::Bool=true,
+                    verbose::Bool=m.options.verbose,
+                    tol::Float64=m.options.tol,
+                    maxiter::Int64=m.options.maxiter,
+                    fctype::FCType=fcgiven,
+                    expectation_horizon::Union{Nothing,Int64}=nothing,
+                    sparse_solver::Function=(A, b) -> A \ b
     )
     NT = length(p.range)
     nauxs = length(m.auxvars)
@@ -113,7 +113,7 @@ function simulate(m::Model, p::Plan, exog_data::AbstractArray{Float64,2};
     end
     if linearize
         org_med = m.evaldata
-        linearize!(m; deviation = deviation)
+        linearize!(m; deviation=deviation)
     end
     if anticipate
         @timer gdata = StackedTimeSolverData(m, p, fctype)
@@ -232,9 +232,9 @@ end
 
 # Simulate command, IRIS style but without a range
 function simulate(m::Model, D1::Dict{<:AbstractString,<:Any}, plan::Plan; 
-    deviation::Bool = false, overlay::Bool = false, kwargs...)::Dict{String,<:Any}
+    deviation::Bool=false, overlay::Bool=false, kwargs...)::Dict{String,<:Any}
     # Convert dictionary to Array{Float64,2}
-    data01 = dict2array(D1, m.varshks, range = plan.range)
+    data01 = dict2array(D1, m.varshks, range=plan.range)
     # Adjust array with steady state values if necessary
     if deviation
         # Check that the steady state has been solved for
@@ -250,11 +250,11 @@ function simulate(m::Model, D1::Dict{<:AbstractString,<:Any}, plan::Plan;
     ig = zeros(0, 0)
     if :initial_guess ∈ keys(kwargs)
         if kwargs[:initial_guess] isa Dict
-            ig = dict2array(kwargs[:initial_guess], m.varshks, range = plan.range)
+            ig = dict2array(kwargs[:initial_guess], m.varshks, range=plan.range)
         end
     end
     # Call native simulate commmand with Array{Float64,2}
-    data02 = simulate(m, plan, data01; kwargs..., initial_guess = ig)
+    data02 = simulate(m, plan, data01; kwargs..., initial_guess=ig)
     # Remove steady state values from data02
     if deviation
         data02 = data02 .- datass;
@@ -271,7 +271,7 @@ function simulate(m::Model, D1::Dict{<:AbstractString,<:Any}, plan::Plan;
 end
 
 # Simulate command, IRIS style with a range
-function simulate(m::Model, D1::Dict{<:AbstractString,<:Any}, rng::AbstractUnitRange, plan::Plan = Plan(m, rng); kwargs...)::Dict{String,<:Any}
+function simulate(m::Model, D1::Dict{<:AbstractString,<:Any}, rng::AbstractUnitRange, plan::Plan=Plan(m, rng); kwargs...)::Dict{String,<:Any}
     # If we have a range, we just take a slice of the plan to enforce the range,
     # but taking into account the model max lag and max lead.
     plan = plan[rng,m];
@@ -281,7 +281,24 @@ function simulate(m::Model, D1::Dict{<:AbstractString,<:Any}, rng::AbstractUnitR
 end
 
 # Simulate command with one date, just in case
-function simulate(m::Model, D1::Dict{<:AbstractString,<:Any}, rng::MIT, plan::Plan = Plan(m, rng:rng); kwargs...)::Dict{String,<:Any}
+function simulate(m::Model, D1::Dict{<:AbstractString,<:Any}, rng::MIT, plan::Plan=Plan(m, rng:rng); kwargs...)::Dict{String,<:Any}
     # return simulate(m,D1,rng:rng, plan; kwargs...)
     return simulate(m, D1, plan[rng:rng,m]; kwargs...)
 end
+
+import ..SimData
+
+simulate(m::Model, D1::SimData, rng, plan::Plan=Plan(m, rng); kwargs...) = simulate(m, D1, plan[rng, m]; kwargs...)
+
+function simulate(m::Model, D1::SimData, plan::Plan;
+    deviation::Bool=false, overlay::Bool=false, kwargs...)::typeof(D1)
+    ret = copy(D1)
+    ig = get(kwargs, :initial_guess, zeros(0, 0))
+    if ig isa SimData
+        ig = ig[plan.range, m.varshks]
+    end
+    sim = simulate(m, plan, D1[plan.range, m.varshks]; kwargs..., initial_guess=ig)
+    # overlay the sim data onto ret
+    ret[plan.range, m.varshks] = sim
+    return ret
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index c070a4e..d55231c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,18 @@
 using StateSpaceEcon
+using TimeSeriesEcon
+using ModelBaseEcon
+
+@info "Loading Model examples. Might take some time to pre-compile."
+@using_example M1
+@using_example M2
+@using_example M3
+@using_example M3nl
+@using_example M6
+@using_example M7
+@using_example M7A
 
 using Test
+using Suppressor
 
 @testset "1dsolvers" begin
     # f(x) = (x-2)*(x-3) = a x^2 + b x + c with vals = [a, x, b, c]
@@ -9,31 +21,50 @@ using Test
         vals = [1.0, NaN, -5.0, 6.0]
 
         vals[2] = 0.0
-        @test StateSpaceEcon.SteadyStateSolver.newton1!(fdf, vals, 2; tol = eps(), maxiter = 8)
+        @test StateSpaceEcon.SteadyStateSolver.newton1!(fdf, vals, 2; tol=eps(), maxiter=8)
         @test vals ≈ [1.0, 2.0, -5.0, 6.0] atol = 1e3 * eps()
 
         vals[2] = 6.0
-        @test StateSpaceEcon.SteadyStateSolver.newton1!(fdf, vals, 2; tol = eps(), maxiter = 8)
+        @test StateSpaceEcon.SteadyStateSolver.newton1!(fdf, vals, 2; tol=eps(), maxiter=8)
         @test vals ≈ [1.0, 3.0, -5.0, 6.0] atol = 1e3 * eps()
 
         vals[2] = 0.0
-        @test StateSpaceEcon.SteadyStateSolver.bisect!(f, vals, 2, fdf(vals)[2][2]; tol = eps())
+        @test StateSpaceEcon.SteadyStateSolver.bisect!(f, vals, 2, fdf(vals)[2][2]; tol=eps())
         @test vals ≈ [1.0, 2.0, -5.0, 6.0] atol = 1e3 * eps()
 
         vals[2] = 6.0
-        @test StateSpaceEcon.SteadyStateSolver.bisect!(f, vals, 2, fdf(vals)[2][2]; tol = eps())
+        @test StateSpaceEcon.SteadyStateSolver.bisect!(f, vals, 2, fdf(vals)[2][2]; tol=eps())
         @test vals ≈ [1.0, 3.0, -5.0, 6.0] atol = 1e3 * eps()
     end
 end
 
-using ModelBaseEcon
-@using_example M1
-@using_example M2
-@using_example M3
-@using_example M3nl
-@using_example M6
-@using_example M7
-@using_example M7A
+@testset "Plans" begin
+    m = M1.model
+    p = Plan(m, 1:3)
+    @test first(p.range) == ii(0)
+    @test last(p.range) == ii(4)
+    out = @capture_out print(p)
+    @test length(split(out, '\n')) == 2
+    @test p[1] == [:y_shk]
+    @test p[ii(1)] == [:y_shk]
+    endogenize!(p, :y_shk, ii(1))
+    @test isempty(p[ii(1)])
+    endogenize!(p, :y_shk, ii(1):ii(3))
+    exogenize!(p, :y, ii(2))
+    exogenize!(p, :y, ii(4))
+    # make sure indexing with integers works as well
+    @test p[ii(0)] == p[1] == [:y_shk]
+    @test p[ii(1)] == p[2] == []
+    @test p[ii(2)] == p[3] == [:y]
+    @test p[ii(3)] == p[4] == []
+    @test p[ii(4)] == p[5] == [:y, :y_shk]
+    out = @capture_out print(p)
+    @test length(split(out, '\n')) == 6
+    out = @capture_out print(IOContext(stdout, :displaysize => (7, 80)), p)
+    @test length(split(out, '\n')) == 4
+    @test length(split(out, '⋮')) == 2
+end
 
+include("simdatatests.jl")
 include("sstests.jl")
 include("simtests.jl")
diff --git a/test/simdatatests.jl b/test/simdatatests.jl
new file mode 100644
index 0000000..acdbde6
--- /dev/null
+++ b/test/simdatatests.jl
@@ -0,0 +1,143 @@
+
+@testset "SimData" begin
+    @test_throws ArgumentError SimData(1M10, (:a, :b, :c), rand(10, 2))
+    let nms = (:a, :b), dta = rand(20, 2), sd = SimData(2000Q1, nms, copy(dta)), dta2 = rand(size(dta)...)
+        @test firstdate(sd) == 2000Q1
+        @test lastdate(sd) == 2000Q1 + 20 - 1
+        @test frequencyof(sd) == Quarterly
+        @test sd isa AbstractMatrix
+        @test size(sd) == size(dta)
+        # integer indexing must be identical to dta
+        for i in axes(dta, 2)
+            @test sd[:,i] == dta[:,i]
+        end
+        for j in axes(dta, 1)
+            @test sd[j,:] == dta[j,:]
+        end
+        for i in eachindex(dta)
+            @test sd[i] == dta[i]
+        end
+        # set indexing with integers
+        for i in eachindex(dta2)
+            sd[i] = dta2[i]
+        end
+        @test sd[:] == dta2[:]
+        for i in axes(dta, 1)
+            sd[i,:] = dta[i,:]
+        end
+        @test sd[:] == dta[:]
+        for i in axes(dta2, 2)
+            sd[:,i] = dta2[:,i]
+        end
+        @test sd[:] == dta2[:]
+        # access by dot notation (for the columns)
+        @test propertynames(sd) == nms
+        @test sd.a isa TSeries
+        @test sd.b isa TSeries
+        @test_throws BoundsError sd.c
+        @test sd.a.values == dta2[:,1]
+        sd.a[:] = dta[:,1]
+        sd.b[:] = dta[:,2]
+        @test sd[:] == dta[:]
+        @test sd[:a] isa TSeries && sd[:a].values == dta[:,1]
+        @test sd["b"] isa TSeries && sd["b"].values == dta[:,2]
+        # 
+        sd.a = dta[:,1]
+        sd.b = dta[:,2]
+        @test sd[:] == dta[:]
+        # 
+        sd[:] = dta[:]
+        sd.a = sd.b
+        @test sd[:,1] == sd[:,2]
+        sd.a = zeros(size(dta, 1))
+        @test sum(abs, sd.a.values) == 0
+        @test_throws DimensionMismatch sd.a = ones(length(sd.a) + 5)
+        # access to rows by MIT
+        sd[:] = dta[:]
+        @test sd[2000Q1] isa NamedTuple{nms,NTuple{length(nms),Float64}}
+        for (i, idx) in enumerate(mitrange(sd))
+            @test [values(sd[idx])...] == dta[i,:]
+        end
+        sd[2000Q1] = zeros(size(dta, 2))
+        @test sum(abs, sd[1,:]) == 0
+        @test_throws BoundsError sd[1999Q4]
+        @test_throws DimensionMismatch sd[2000Q3] = zeros(size(dta, 2) + 1)
+        sd[2000Q2] = (b = 7.0, a = 1.5)
+        @test sd[2,:] == [1.5, 7.0]
+        sd[2004Q4] = (a = 5.0,)
+        @test sd[end,1] == 5.0
+        @test_throws BoundsError sd[2004Q4] = (a = 5.0, c = 1.1)
+        @test_throws ArgumentError sd[2000Y]
+        @test_throws BoundsError sd[1999Q4] = (a = 5.0, b = 1.1)
+        @test_throws BoundsError sd[2010Q4] = rand(2)
+        # 
+        @test sd[2001Q1:2002Q4] isa SimData
+        @test sd[2001Q1:2002Q4] == sd[5:12,:]
+        sd[2001Q1:2002Q4] = 1:16
+        @test sd[5:12,:] == reshape(1.0:16.0, :, 2)
+        @test_throws BoundsError sd[1111Q4:2002Q1]
+        @test_throws BoundsError sd[2002Q1:2200Q2]
+        @test_throws DimensionMismatch sd[2001Q1:2002Q4] = 1:17
+        # assign new column
+        sd1 = hcat(sd, c=sd.a + 3.0)
+        @test sd1[nms] == sd
+        @test sd1[[:a,:c]] == sd1[:,[1,3]]
+        # access with 2 args MIT and Symbol
+        @test sd[2001Q2, (:a, :b)] isa NamedTuple
+        let foo = sd[2001Q2:2002Q1, (:a, :b)] 
+            @test foo isa SimData
+            @test size(foo) == (4, 2) 
+            @test firstdate(foo) == 2001Q2
+        end
+        @test_throws BoundsError sd[1999Q1, (:a,)]
+        @test_throws BoundsError sd[2001Q1:2001Q2, (:a, :c)]
+        @test_throws Exception sd.c = 5
+        @test similar(sd) isa typeof(sd)
+        @test_throws BoundsError sd[1999Q1:2000Q4] = zeros(8, 2)
+        @test_throws BoundsError sd[2004Q1:2013Q4, [:a, :b]]
+        # setindex with two 
+        sd[2001Q1, (:b,)] = 3.5
+        @test sd[5,2] == 3.5
+        sd[2000Q1:2001Q4, (:b,)] = 3.7
+        @test all(sd[1:8,2] .== 3.7)
+        @test_throws BoundsError sd[1999Q1:2000Q4, (:a, :b)] = 5.7
+        @test_throws BoundsError sd[2000Q1:2000Q4, (:a, :c)] = 5.7
+    end
+end
+
+
+@testset "SimData show" begin
+    # test the case when column labels are longer than the numbers
+    let io = IOBuffer(), sd =  SimData(1U, (:verylongandsuperboringnameitellya, :anothersuperlongnamethisisridiculous, :a), rand(20, 3) .* 100)
+        show(io, sd)
+        lines = readlines(seek(io, 0))
+        # labels longer than 10 character are abbreviated with '…' at the end
+        @test length(split(lines[2], '…')) == 3
+    end
+    let io = IOBuffer() , sd = SimData(1U, (:alpha, :beta), zeros(24, 2))
+        show(io, sd)
+        lines = readlines(seek(io, 0))
+        lens = length.(lines)
+        # when labels are longer than the numbers, the entire column stretches to fit the label
+        @test lens[2] == lens[3]
+    end
+    nrow = 24
+    letters = Symbol.(['a':'z'...])
+    for (nlines, fd) = zip([3, 4, 5, 6, 7, 8, 22, 23, 24, 25, 26, 30], Iterators.cycle((qq(2010, 1), mm(2010, 1), yy(2010), ii(1))))
+        for ncol = [2,5,10,20]
+            # display size if nlines × 80
+            # data size is nrow × ncol
+            # when printing data there are two header lines - summary and column names
+            io = IOBuffer()
+            sd = SimData(fd, tuple(letters[1:ncol]...), rand(nrow, ncol))
+            show(IOContext(io, :displaysize => (nlines, 80)), MIME"text/plain"(), sd)
+            lines = readlines(seek(io, 0))
+            @test length(lines) == max(3, min(nrow + 2, nlines - 3))
+            @test maximum(length, lines) <= 80
+            io = IOBuffer()
+            show(IOContext(io, :limit => false), sd)
+            lines = readlines(seek(io, 0))
+            @test length(lines) == nrow + 2
+        end
+    end
+end
\ No newline at end of file
diff --git a/test/simtests.jl b/test/simtests.jl
index 914cea5..3f01596 100644
--- a/test/simtests.jl
+++ b/test/simtests.jl
@@ -7,10 +7,10 @@
     m.β = 1.0 - m.α
     for T = 6:10
         p = Plan(m, 2:T - 1)
-        data = zeroarray(m, p)
+        data = zerodata(m, p)
         data[1, :] = [1 0]   # initial condition
         data[end, :] = [5 0] # final condition
-        sim00 = simulate(m, p, data)
+        sim00 = simulate(m, data, p)
         exp00 = hcat(1.0:(5.0 - 1.0) / (T - 1):5.0, zeros(T))
         # @info "T=$T" sim00 exp00 data
         @test sim00 ≈ exp00
@@ -18,18 +18,18 @@
         data1 = copy(data)
         shk = .1
         y2val = ((T - 2) * (1 + 2 * shk) + 5.0) / (T - 1)
-        data1[2,2] = shk  # shock at time 2
-        sim01 = simulate(m, p, data1)
+        data1[2U,:y_shk] = shk  # shock at time 2
+        sim01 = simulate(m, data1, p)
         exp01 = vcat([1.0 y2val:(5.0 - y2val) / (T - 2):5.0...], [0 shk zeros(T - 2)...])'
         @test sim01 ≈ exp01
         # exogenous-endogenous swap 
         # - replicate the solution above backing out the shock that produces it
         data2 = copy(data)
-        data2[2, 1] = sim01[2, 1]
+        data2[2U, :y] = sim01[2U, :y]
         p2 = deepcopy(p)
         exogenize!(p2, :y, 2)
         endogenize!(p2, :y_shk, 2)
-        sim02 = simulate(m, p2, data2)
+        sim02 = simulate(m, data2, p2)
         exp02 = sim01
         @test sim02 ≈ exp02
     end
@@ -139,8 +139,8 @@ end
     data[1,1] = 1.0;
     data[end,1] = 5.0;
     data[3,2] = 0.1;
-    res_u = simulate(m, p, data; anticipate = false)
-    true_u = Float64[1 0; 2 0; 47/15 0.1; 61/15 0; 5 0]
+    res_u = simulate(m, p, data; anticipate=false)
+    true_u = Float64[1 0; 2 0; 47 / 15 0.1; 61 / 15 0; 5 0]
     @test res_u ≈ true_u atol = 1e-12
 end
 
@@ -171,7 +171,7 @@ end
             # Run the simulations and test
             ares = simulate(m, p, adata)
             @test ares ≈ atrue atol = m.options.tol
-            ures = simulate(m, p, udata; anticipate = false)
+            ures = simulate(m, p, udata; anticipate=false)
             @test ures ≈ utrue atol = m.options.tol
         end
         let adata = zeroarray(m, p),
@@ -184,7 +184,7 @@ end
             # Run the simulations and test
             ares = simulate(m, p, adata)
             @test ares ≈ atrue atol = m.options.tol
-            ures = simulate(m, p, udata; anticipate = false)
+            ures = simulate(m, p, udata; anticipate=false)
             @test ures ≈ utrue atol = m.options.tol
         end
     end
@@ -210,7 +210,7 @@ end
         # Run the simulations and test
         ares = simulate(m, p, adata)
         @test ares ≈ atrue atol = m.options.tol
-        ures = simulate(m, p, udata; anticipate = false)
+        ures = simulate(m, p, udata; anticipate=false)
         @test ures ≈ utrue atol = m.options.tol
     end
 end
@@ -242,7 +242,7 @@ end
             # Run the simulations and test
             ares = simulate(m, p, adata)
             @test ares ≈ atrue atol = m.options.tol
-            ures = simulate(m, p, udata; anticipate = false)
+            ures = simulate(m, p, udata; anticipate=false)
             @test ures ≈ utrue atol = m.options.tol
         end
         let adata = zeroarray(m, p),
@@ -255,7 +255,7 @@ end
             # Run the simulations and test
             ares = simulate(m, p, adata)
             @test ares ≈ atrue atol = m.options.tol
-            ures = simulate(m, p, udata; anticipate = false)
+            ures = simulate(m, p, udata; anticipate=false)
             @test ures ≈ utrue atol = m.options.tol
         end
     end
@@ -281,7 +281,7 @@ end
         # Run the simulations and test
         ares = simulate(m, p, adata)
         @test ares ≈ atrue atol = m.options.tol
-        ures = simulate(m, p, udata; anticipate = false)
+        ures = simulate(m, p, udata; anticipate=false)
         @test ures ≈ utrue atol = m.options.tol
     end
 end