diff --git a/README.md b/README.md
index 75a43d241..72d1292a3 100644
--- a/README.md
+++ b/README.md
@@ -255,7 +255,9 @@ Configure and make:
     - Next to `Command Line Tools`, click on `Install`
 - Download [MacOSX10.6.sdk.zip](http://www.jamesgeorge.org/uploads/MacOSX10.6.sdk.zip) and unzip it in `/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs`. Change the owner to be `root` and the group to be `wheel`
 
-- Install [Homebrew](http://mxcl.github.io/homebrew/) and [Homebrew Science](https://github.com/Homebrew/homebrew-science)
+- Install [Homebrew](http://mxcl.github.io/homebrew/) by following the instructions on the website
+- Tap [Homebrew Science](https://github.com/Homebrew/homebrew-science) by doing:
+    - ```brew tap homebrew/science```
 - Install the following brews:
     - ```brew install automake```
     - ```brew install gsl```
diff --git a/doc/dseries-and-reporting/dseriesReporting.tex b/doc/dseries-and-reporting/dseriesReporting.tex
index 441888cca..24b2ef8b2 100644
--- a/doc/dseries-and-reporting/dseriesReporting.tex
+++ b/doc/dseries-and-reporting/dseriesReporting.tex
@@ -104,7 +104,7 @@ X =
     \texttt{multiplyByTwo()}.
     \myitem In most object-oriented
     languages, writing \texttt{X.multiplyByTwo();} will change the value
-    contained in \texttt{X} to \texttt{2/2}
+    contained in \texttt{X} to \texttt{2}
     \end{itemize}
   \end{itemize}
 \end{frame}
@@ -139,7 +139,7 @@ X =
     \item[] \begin{itemize}
     \myitem To get the desired change, you must overwrite \texttt{X}
 \begin{verbatim}
->> X=X.multiplyByTwo()
+>> X = X.multiplyByTwo()
 
 X =
 
@@ -194,10 +194,6 @@ X =
   >> t = dates(`1999q1', `2020q2', `-190q3');
 \end{alltt}
 Notice that noncontiguous and negative dates are possible
-    \myitem{A date range}
-\begin{alltt}
-  >> dr = dates(`1999y'):dates(`2020y');
-\end{alltt}
     \myitem Can also create \texttt{dates} programatically
 \begin{alltt}
   >> t = dates(4, [1990; 1990; 1978], [1; 2; 3])
@@ -234,10 +230,6 @@ Notice that noncontiguous and negative dates are possible
     \myitem{Multiple dates} \texttt{t = [1999q1 2020q2 1960q3];}
 \begin{alltt}
   \(\Rightarrow\) t = [dates(`1999q1') dates(`2020q2') dates(`1960q3')];
-\end{alltt}
-    \myitem{A date range} \texttt{dr = 1999y:2020y;}
-\begin{alltt}
-  \(\Rightarrow\) dr = dates(`1999y'):dates(`2020y');
 \end{alltt}
     \myitem NB: This can cause problems when dates are included in strings. \textit{e.g.,}
 \begin{alltt}
@@ -398,7 +390,7 @@ would be transformed into
     \myitem All time series in a \texttt{dseries} must have the same frequency
     \myitem A \texttt{dseries} runs from the earliest date to the latest date,
     with \texttt{NaN}'s inserted to pad the shorter series
-    \myitem A \texttt{dseries} object contains 3 members:
+    \myitem A \texttt{dseries} object contains 6 members:
     \begin{itemize}
       \myitem{\textbf{\texttt{data}}}: The data points
       \myitem{\textbf{\texttt{dates}}}: The \texttt{dates} of the sample
@@ -452,7 +444,7 @@ would be transformed into
     \texttt{NAMES\_\_}, and, optionally, \texttt{TEX\_\_}. More info in the
     manual. Data are vectors.
     \begin{alltt}
-      INIT__ = `1999q1’;
+      INIT__ = `1999q1';
       NAMES__ = \{`cons'\};
       cons = randn(100,1);
     \end{alltt}
@@ -475,13 +467,11 @@ would be transformed into
     \begin{alltt}
   ts(2001Y:2003Y)
     \end{alltt}
-    \myitem Can also use integer indices
+    \myitem Can also use integer indices (in a roundabout way)
     \begin{alltt}
   start = find(ts.dates==2001Y);
-  ts(start:end);
+  ts(ts.dates(start:end));
     \end{alltt}
-    \myitem NB: Integer indices cannot be used to extract a single observation
-    because of lag/lead notation (see below)
   \end{itemize}
 \end{frame}
 
@@ -687,6 +677,12 @@ the time range of \verb+vs+ will be the union of \verb+ts.dates+ and
 \end{frame}
 
 
+\tikzstyle{abstract}=[rectangle, draw=black, rounded corners, fill=blue!10, drop shadow,
+        text centered, anchor=north, text=black, text width=2.6cm]
+\tikzstyle{comment}=[rectangle, draw=black, rounded corners, fill=green!10, drop shadow,
+        text centered, anchor=north, text=black, text width=2.6cm]
+
+
 \begin{frame}[fragile,t]
   \frametitle{Reporting Class Hierarchy}
   \begin{itemize}
@@ -694,9 +690,9 @@ the time range of \verb+vs+ will be the union of \verb+ts.dates+ and
   \myitem Arrows represent what the new object can be added to; objects in ellipses are treated a bit differently (explained below)
   \end{itemize}
   \begin{center}
-    \scriptsize{
+    \footnotesize{
     \begin{tikzpicture}[
-        node distance = .45cm,
+        node distance = .25cm,
         auto,
         line/.style={->, >=stealth'},
       ]
@@ -715,29 +711,29 @@ the time range of \verb+vs+ will be the union of \verb+ts.dates+ and
               \textbf{Section}
               \nodepart{second}\texttt{addSection(...);}
             };
-      \node (Paragraph) [abstract, ellipse split, rectangle split parts=2, right=of Section]
-            {
-              \textbf{Paragraph}
-              \nodepart{lower}\texttt{addParagraph(...);}
-            };
-      \node (Vspace) [abstract, ellipse split, rectangle split parts=2, left=of Section]
-            {
-              \textbf{Vspace}
-              \nodepart{lower}\texttt{addVspace(...);}
-            };
-      \node (Spacer) [abstract, rectangle split, rectangle split parts=2, opacity=0, below=of Section]
+      \node (Spacer) [rectangle split, rectangle split parts=2, opacity=0, below=of Section]
             {
             };
-      \node (Graph) [abstract, rectangle split, rectangle split parts=2, left=of Spacer]
+      \node (Graph) [abstract, rectangle split, rectangle split parts=2, left=of Spacer, xshift=.25cm]
             {
               \textbf{Graph}
               \nodepart{second}\texttt{addGraph(...);}
             };
-      \node (Table) [abstract, rectangle split, rectangle split parts=2, right=of Spacer, text height=]
+      \node (Table) [abstract, rectangle split, rectangle split parts=2, right=of Spacer, xshift=-.25cm, text height=]
             {
               \textbf{Table}
               \nodepart{second}\texttt{addTable(...);}
             };
+      \node (Vspace) [comment, rectangle split, rectangle split parts=2, left=of Graph]
+            {
+              \textbf{Vspace}
+              \nodepart{second}\texttt{addVspace(...);}
+            };
+      \node (Paragraph) [comment, rectangle split, rectangle split parts=2, right=of Table]
+            {
+              \textbf{Paragraph}
+              \nodepart{second}\texttt{addParagraph(...);}
+            };
       \node (Series) [abstract, rectangle split, rectangle split parts=2, below=of Spacer]
             {
               \textbf{Series}
@@ -747,8 +743,8 @@ the time range of \verb+vs+ will be the union of \verb+ts.dates+ and
       \draw [line] (Series) to node { } (Graph);
       \draw [line] (Table) to node { } (Section);
       \draw [line] (Graph) to node { } (Section);
-      \draw [line] (Paragraph) to node { } (Section);
-      \draw [line] (Vspace) to node { } (Section);
+      \draw [line] (Paragraph.north) to node { } (Section);
+      \draw [line] (Vspace.north) to node { } (Section);
       \draw [line] (Section) to node { } (Page);
       \draw [line] (Page) to node { } (Report);
     \end{tikzpicture}}
@@ -819,8 +815,8 @@ the time range of \verb+vs+ will be the union of \verb+ts.dates+ and
       \myitem \textbf{Options}: \texttt{data}, \texttt{graphDirName}, \texttt{graphName}, \texttt{graphSize}, \texttt{height}, \ldots
       \begin{alltt}
 >> rep = rep.addGraph(`title', `Headline Inflation (y/y)', ...
-                      `xrange', prange, ...
-                      `shade', srange, ...
+                      `xrange', dates(`2007q1'):dates(`2013q4'), ...
+                      `shade', dates(`2010q1'):dates(`2013q4'), ...
                       `showZeroline', true);
       \end{alltt}
     \end{itemize}
@@ -829,7 +825,7 @@ the time range of \verb+vs+ will be the union of \verb+ts.dates+ and
       \myitem \textbf{Options}: \texttt{data}, \texttt{showHlines}, \texttt{precision}, \texttt{range}, \texttt{seriesToUse}, \ldots
       \begin{alltt}
 >> rep = rep.addTable(`title', \{`Real GDP Growth',`subtitle 1'\}, ...
-                      `range', larange, ...
+                      `range', dates(`2007y'):dates(`2014y'), ...
                       `vlineAfter', dates(`2011y'));
       \end{alltt}
     \end{itemize}
@@ -852,7 +848,7 @@ the time range of \verb+vs+ will be the union of \verb+ts.dates+ and
   To create a report:
   \begin{itemize}
   \myitem \texttt{write()}: Writes the report to a \LaTeX\ file
-  \myitem \texttt{compile(\ldots)}: Compiles the report
+  \myitem \texttt{compile(\ldots)}: Compiles the report, creating a \texttt{.pdf} file
     \begin{itemize}
     \myitem \textbf{Options}: \texttt{compiler}
     \end{itemize}
@@ -867,43 +863,49 @@ the time range of \verb+vs+ will be the union of \verb+ts.dates+ and
 \end{frame}
 
 
-\subsection{Examples}
-
 \section{Putting it All Together}
 \begin{frame}[fragile=singleslide,t]
-  \frametitle{Create Report of IRFs from \texttt{example1.mod}}
+  \frametitle{Create Report of IRFs from \texttt{example1.mod} (1/3)}
   \begin{itemize}
   \myitem \texttt{example1.mod} is located in the Dynare \texttt{examples} directory
   \myitem The lines below can be added at the end of that file.
   \end{itemize}
   \begin{block}{Create \texttt{dseries} from IRFs}
 \begin{verbatim}
+shocke = dseries();
+shocku = dseries();
 @#define endovars=["y", "c", "k", "a", "h", "b"]
 @#for var in endovars
-  shocke.@{var} = dseries(@{var}_e, 2014q3, `@{var}');
-  shocku.@{var} = dseries(@{var}_u, 2014q3, `@{var}');
+  shocke = [shocke dseries(@{var}_e, 2014q3, '@{var}')];
+  shocku = [shocku dseries(@{var}_u, 2014q3, '@{var}')];
 @#endfor
 \end{verbatim}
   \end{block}
 \end{frame}
 
 \begin{frame}[fragile=singleslide,t]
-  \frametitle{Create Report of IRFs from \texttt{example1.mod}}
-  \begin{block}{Populate Report (1/2)}
+  \frametitle{Create Report of IRFs from \texttt{example1.mod} (2/3)}
+  \begin{block}{Populate Report}
 \small{
 \begin{verbatim}
+r = report();
 @#for shock in ["e", "u"]
-    report = report.addPage(`title', {`Dseries \& Report Example', ...
-                                      `Shock to @{shock}'}, ...
-                            `titleFormat', {`\Large\bfseries', ...
-                                            `\large\bfseries'});
-    report = report.addSection(`cols', 2);
+    r = r.addPage(`title',{`Dseries/Report Example',`Shock @{shock}'},...
+                  `titleFormat', {`\Large\bfseries', `\large\bfseries'});
+    r = r.addSection(`cols', 2);
 @#  for var in endovars
-      report = report.addGraph(`data', shock@{shock}.@{var}, ...
-                               `title', `@{var}', ...
-                               `showGrid', false, ...
-                               `showZeroLine', true);
+      r = r.addGraph(`data', shock@{shock}.@{var}, `title', `@{var}', ...
+                     `showGrid', false, `yTickLabelPrecision', 2, ...
+                     `yTickLabelZeroFill', false);
+      r = r.addSeries(`graphHline', 0, `graphLineColor', `red');
 @#  endfor
+    r = r.addVspace(`number', 2);
+    r = r.addSection(`cols', 1);
+    r = r.addTable(`range', 2022q1:2024q1, `precision', 5);
+@#  for var in endovars
+      r = r.addSeries(`data', shock@{shock}.@{var});
+@#  endfor
+@#endfor
 \end{verbatim}
 }
   \end{block}
@@ -911,27 +913,29 @@ the time range of \verb+vs+ will be the union of \verb+ts.dates+ and
 
 
 \begin{frame}[fragile=singleslide,t]
-  \frametitle{Create Report of IRFs from \texttt{example1.mod}}
-  \begin{block}{Populate Report (2/2)}
-\small{
-\begin{verbatim}
-    report = report.addVspace(`number', 2);
-    report = report.addSection(`cols', 1);
-    report = report.addTable(`range', 2022q1:2024q1, `precision', 5);
-
-@#  for var in endovars
-      report = report.addSeries(`data', shock@{shock}.@{var});
-@#  endfor
-@#endfor
-\end{verbatim}
-}
-  \end{block}
-
+  \frametitle{Create Report of IRFs from \texttt{example1.mod} (3/3)}
   \begin{block}{Compile Report}
 \small{
 \begin{verbatim}
-report.write();
-report.compile();
+r.write();
+r.compile();
+\end{verbatim}
+}
+  \end{block}
+  \begin{block}{Output Files}
+\small{
+\begin{verbatim}
+>> ls report.*
+report.aux  report.log  report.pdf  report.synctex.gz  report.tex
+
+>> ls tmpRepDir/
+graph_pg1_sec1_row1_col1.tex  graph_pg2_sec1_row1_col1.tex
+graph_pg1_sec1_row1_col2.tex  graph_pg2_sec1_row1_col2.tex
+graph_pg1_sec1_row2_col1.tex  graph_pg2_sec1_row2_col1.tex
+graph_pg1_sec1_row2_col2.tex  graph_pg2_sec1_row2_col2.tex
+graph_pg1_sec1_row3_col1.tex  graph_pg2_sec1_row3_col1.tex
+graph_pg1_sec1_row3_col2.tex  graph_pg2_sec1_row3_col2.tex
+table_pg1_sec2_row1_col1.tex  table_pg2_sec2_row1_col1.tex
 \end{verbatim}
 }
   \end{block}
diff --git a/doc/dynare.texi b/doc/dynare.texi
index da47c7a88..d0aed6a71 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -9713,12 +9713,6 @@ A @code{nobs}*1 cell of strings or a @code{nobs}*p character array, the names of
 @item tex
  A @code{nobs}*1 cell of strings or a @code{nobs}*p character array, the tex names of the variables.
 
-@item freq
-A scalar integer equal to 1, 4, 12 or 52, the frequency of the dataset.
-
-@item init
-A single element @dates object, the initial date of the sample.
-
 @item dates
 A @dates object with @code{nobs} element, the dates of the sample.
 
@@ -9727,13 +9721,13 @@ A @code{nobs} by @code{vobs} array of doubles, the data.
 
 @end table
 
-@noindent @code{freq}, @code{nobs}, @code{vobs}, @code{data}, @code{name}, @code{tex} are private members. The following constructors are available:
+@noindent @code{nobs}, @code{vobs}, @code{data}, @code{name}, @code{tex} are private members. The following constructors are available:
 
 @deftypefn  {dseries} dseries ()
 @deftypefnx {dseries} dseries (@var{INITIAL_DATE})
 @deftypefnx {dseries} dseries (@var{RANGE_OF_DATES})
 
-Instantiates an empty @dseries object, with, if defined, an initial date given by the single element @dates object @var{INITIAL_DATE} or the first element of the @dates object @var{RANGE_OF_DATES} (the frequency is then set accordingly).
+Instantiates an empty @dseries object, with, if defined, an initial date given by the single element @dates object @var{INITIAL_DATE} or the first element of the @dates object @var{RANGE_OF_DATES}.
 
 @end deftypefn
 
@@ -10160,6 +10154,24 @@ ans is a dseries object:
 
 @sp 1
 
+@deftypefn {dseries} {@var{f} =} freq (@var{B})
+
+Returns the frequency of the variables in @dseries object @var{B}.
+
+@examplehead
+@example
+>> ts = dseries(randn(3,2),'1973Q1');
+>> ts.freq
+
+ans =
+
+     4
+@end example
+
+@end deftypefn
+
+@sp 1
+
 @deftypefn{dseries} {@var{D} =} horzcat (@var{A}, @var{B}[, ...])
 
 Overloads the @code{horzcat} Matlab/Octave's method for @dseries
@@ -10277,6 +10289,21 @@ The previous code should produce something like:
 
 @sp 1
 
+@deftypefn {dseries} {@var{f} =} init (@var{B})
+
+Returns the initial date in @dseries object @var{B}.
+
+@examplehead
+@example
+>> ts = dseries(randn(3,2),'1973Q1');
+>> ts.init
+ans = <dates: 1973Q1>
+@end example
+
+@end deftypefn
+
+@sp 1
+
 @deftypefn {dseries} {@var{C} = } insert (@var{A}, @var{B}, @var{I})
 
 Inserts variables contained in @dseries object @var{B} in @dseries object @var{A} at positions specified by integer scalars in vector @var{I}, returns augmented @dseries object @var{C}. The integer scalars in @var{I} must take values between @code{1} and @code{A.length()+1} and refers to @var{A}'s column numbers. The @dseries objects @var{A} and @var{B} need not to be defined over the same time ranges, but it is assumed that they have common frequency.
@@ -10393,6 +10420,21 @@ ans is a dseries object:
 
 @sp 1
 
+@deftypefn {dseries} {@var{l} =} last (@var{B})
+
+Returns the last date in @dseries object @var{B}.
+
+@examplehead
+@example
+>> ts = dseries(randn(3,2),'1973Q1');
+>> ts.last
+ans = <dates: 1973Q3>
+@end example
+
+@end deftypefn
+
+@sp 1
+
 @deftypefn {dseries} {@var{B} = } lead (@var{A}[, @var{p}])
 
 Returns leaded time series. Default value of @var{p}, the number of leads, is @code{1}. As for the @code{lag} method, the @dseries class overloads the parenthesis so that @code{ts.lead(p)} is equivalent to @code{ts(p)}.
@@ -10848,6 +10890,40 @@ ts1 is a dseries object:
 
 @sp 1
 
+@deftypefn{dseries} {@var{C} =} remove (@var{A}, @var{B})
+
+Alias for the @code{pop} method with two arguments. Removes variable @var{B} from @dseries object @var{A}.
+
+@examplehead
+@example
+>> ts0 = dseries(ones(3,3));
+>> ts1 = ts0.remove('Variable_2');
+
+ts1 is a dseries object:
+
+   | Variable_1 | Variable_3
+1Y | 1          | 1
+2Y | 1          | 1
+3Y | 1          | 1
+@end example
+
+@sp 1
+
+A shorter syntax is available: @code{remove(ts,'Variable_2')} is
+equivalent to @code{ts@{'Variable_2'@} = []} (@code{[]} can be replaced
+by any empty object). This alternative syntax is usefull if more than
+one variable has to be removed. For instance:
+@example
+ts@{'Variable_@@2,3,4@@'@} = [];
+@end example
+will remove @code{Variable_2}, @code{Variable_3} and @code{Variable_4}
+from @dseries object @code{ts} (if these variables exist). Regular
+expressions cannot be used but implicit loops can.
+
+@end deftypefn
+
+@sp 1
+
 @deftypefn{dseries} {@var{B} =} rename (@var{A},@var{oldname},@var{newname})
 
 Rename variable @var{oldname} to @var{newname} in @dseries object
@@ -11163,7 +11239,7 @@ command. Default: @code{`!'}
 @end defmethod
 
 @anchor{addGraph}
-@defmethod Report addGraph data, graphDirName, graphName, graphSize, height, showGrid, showLegend, showLegendBox, legendLocation, legendOrientation, legendFontSize, miscTikzAxisOptions, miscTikzPictureOptions, seriesToUse, shade, shadeColor, shadeOpacity, title, width, xlabel, ylabel, xAxisTight, xrange, xTicks, xTickLabels, xTickLabelAnchor, xTickLabelRotation, yAxisTight, yrange, showZeroline
+@defmethod Report addGraph data, graphDirName, graphName, graphSize, height, showGrid, showLegend, showLegendBox, legendLocation, legendOrientation, legendFontSize, miscTikzAxisOptions, miscTikzPictureOptions, seriesToUse, shade, shadeColor, shadeOpacity, title, width, xlabel, ylabel, xAxisTight, xrange, xTicks, xTickLabels, xTickLabelAnchor, xTickLabelRotation, yAxisTight, yTickLabelFixed, yTickLabelPrecision, yTickLabelZeroFill, yrange, showZeroline
 Adds a @code{Graph} to a @code{Section}.
 @optionshead
 @table @code
@@ -11303,6 +11379,17 @@ The boundary on the y-axis to display in the graph, represented as a
 @code{NUMERICAL_VECTOR} of size @math{2}, with the first entry less
 than the second entry. Default: all
 
+@item yTickLabelFixed, @code{BOOLEAN}
+Round the y tick labels to a fixed number of decimal places, given by
+@ref{yTickLabelPrecision}. Default: @code{true}
+
+@anchor{yTickLabelPrecision}
+@item yTickLabelPrecision, @code{INTEGER}
+The precision with which to report the yTickLabel. Default: @code{1}
+
+@item yTickLabelZeroFill, @code{BOOLEAN}
+Whether or not to fill missing precision spots with zeros. Default: @code{true}
+
 @item showZeroline, @code{BOOLEAN}
 Display a solid black line at @math{y = 0}. Default: @code{false}
 
diff --git a/matlab/@dates/colon.m b/matlab/@dates/colon.m
index 4449d9a3f..a68b82519 100644
--- a/matlab/@dates/colon.m
+++ b/matlab/@dates/colon.m
@@ -30,7 +30,7 @@ function C = colon(varargin) % --*-- Unitary tests --*--
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-if isequal(nargin,2)    
+if isequal(nargin,2)
     A = varargin{1};
     B = varargin{2};
     d = 1;
@@ -68,22 +68,45 @@ end
 
 C = dates();
 n = (B-A)+1;
-C.freq = A.freq;
-C.ndat = n;
-C.time = NaN(n,2);
-C.time(1,:) = A.time;
-
-current_date = A;
-linee = 1;
-
-while current_date<B
-    linee = linee+1;
-    C.time(linee,:) = add_periods_to_array_of_dates(C.time(linee-1,:), C.freq, d);
-    current_date = current_date + d;
+m = n;
+if d>1
+    m = length(1:d:n);
 end
+C.freq = A.freq;
 
-C.time = C.time(1:linee,:);
-C.ndat = rows(C.time);
+if isequal(C.freq,1)
+    C.ndat = m;
+    C.time = NaN(m,2);
+    C.time(:,1) = A.time(1)+transpose(0:d:n-1);
+    C.time(:,2) = 1;
+else
+    C.time = NaN(n,2);
+    initperiods = min(C.freq-A.time(2)+1,n);
+    C.time(1:initperiods,1) = A.time(1);
+    C.time(1:initperiods,2) = transpose(A.time(2)-1+(1:initperiods));
+    if n>initperiods
+        p = n-initperiods;
+        if p<=C.freq
+            C.time(initperiods+(1:p),1) = A.time(1)+1;
+            C.time(initperiods+(1:p),2) = transpose(1:p);
+        else
+            q = fix(p/C.freq);
+            r = rem(p,C.freq);
+            C.time(initperiods+(1:C.freq*q),2) = repmat(transpose(1:C.freq),q,1);
+            C.time(initperiods+(1:C.freq*q),1) = kron(A.time(1)+transpose(1:q),ones(C.freq,1));
+            if r>0
+                C.time(initperiods+C.freq*q+(1:r),1) = C.time(initperiods+C.freq*q,1)+1;
+                C.time(initperiods+C.freq*q+(1:r),2) = transpose(1:r);
+            end
+        end
+    end
+    if d>1
+        C.time = C.time(1:d:n,:);
+        C.ndat = m;
+    else
+        C.ndat = n;
+    end
+end
 
 %@test:1
 %$ % Define two dates
@@ -149,4 +172,21 @@ C.ndat = rows(C.time);
 %$ t(1) = dyn_assert(d.time,e.time);
 %$ t(2) = dyn_assert(d.freq,e.freq);
 %$ T = all(t);
-%$ @eof:3
\ No newline at end of file
+%$ @eof:3
+
+%$ @test:4
+%$ % Create an empty dates object for quaterly data
+%$ qq = dates('Q');
+%$
+%$ % Define expected results.
+%$ e.freq = 4;
+%$ e.time = [1950 1; 1950 2; 1950 3];
+%$
+%$ % Call the tested routine.
+%$ d = qq(1950,1):qq(1950,3);
+%$
+%$ % Check the results.
+%$ t(1) = dyn_assert(d.time,e.time);
+%$ t(2) = dyn_assert(d.freq,e.freq);
+%$ T = all(t);
+%$ @eof:4
\ No newline at end of file
diff --git a/matlab/@dates/dates.m b/matlab/@dates/dates.m
index 8ba7efd3b..99093bfe7 100644
--- a/matlab/@dates/dates.m
+++ b/matlab/@dates/dates.m
@@ -60,10 +60,16 @@ function dd = dates(varargin) % --*-- Unitary tests --*--
 % GNU General Public License for more details.
 %
 % You should have received a copy of the GNU General Public License
- % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-dd = struct('ndat', 0, 'freq', NaN(0), 'time', NaN(0,2));
-dd = class(dd,'dates');
+if nargin>0 && ischar(varargin{1}) && isequal(varargin{1},'initialize')
+    dd = struct('ndat', 0, 'freq', NaN(0), 'time', NaN(0,2));
+    dd = class(dd,'dates');
+    assignin('base','emptydatesobject',dd);
+    return
+end
+
+dd = evalin('base','emptydatesobject');
 
 switch nargin
   case 0
@@ -104,7 +110,7 @@ switch nargin
         dd = subsref(dd,S);
         return
     else
-        error(['dates::dates: Wrong calling sequence!'])
+        error('dates::dates: Wrong calling sequence!')
     end
     for i=2:dd.ndat
         if isdate(varargin{i})
@@ -112,10 +118,10 @@ switch nargin
             if isequal(date.freq,dd.freq)
                 dd.time(i,:) = date.time;
             else
-                 error(['dates::dates: Check that all the inputs have the same frequency (see input number ' str2num(i) ')!'])
+                 error(sprintf('dates::dates: Check that all the inputs have the same frequency (see input number %i)!',i))
             end
         else
-            error(['dates::dates: Input ' str2num(i) ' has to be a string date!'])
+            error(sprintf('dates::dates: Input %i has to be a string date!',i))
         end
     end
 end
diff --git a/matlab/@dates/ge.m b/matlab/@dates/ge.m
index 8f18cf7ab..625f1bbdf 100644
--- a/matlab/@dates/ge.m
+++ b/matlab/@dates/ge.m
@@ -43,26 +43,41 @@ if isequal(A.ndat, B.ndat)
     C = (A==B);
     idx = find(C==0);
     for i=1:length(idx)
-        C(idx(i)) = compare_vectors(@gt, A.time(idx(i),:), B.time(idx(i),:));
+        C(idx(i)) = greaterorequal(A.time(idx(i),:), B.time(idx(i),:));
     end
 else
-    if isequal(A.ndat,1) && isequal(B.ndat,1)
-        C = compare_vectors(@ge, A.time, B.time);
-    elseif isequal(A.ndat,1)
+    if isequal(A.ndat,1)
         C = NaN(B.ndat,1);
         for i=1:B.ndat
-            C(i) = compare_vectors(@ge, A.time, B.time(i,:));
+            C(i) = greaterorequal(A.time, B.time(i,:));
         end
     elseif isequal(B.ndat,1)
         C = NaN(A.ndat,1);
         for i=1:A.ndat
-            C(i) = compare_vectors(@ge, A.time(i,:), B.time);
+            C(i) = greaterorequal(A.time(i,:), B.time);
         end
     else
         C = 0;
     end
 end
 
+
+function c = greaterorequal(a,b)
+    if a(1)>b(1)
+        c = 1;
+    else
+        if a(1)<b(1)
+            c = 0;
+        else
+            if a(2)>=b(2)
+                c = 1;
+            else
+                c = 0;
+            end
+        end
+    end
+
+
 %@test:1
 %$ % Define some dates
 %$ date_2 = '1950Q2';
diff --git a/matlab/@dates/gt.m b/matlab/@dates/gt.m
index 241dd729e..796038885 100644
--- a/matlab/@dates/gt.m
+++ b/matlab/@dates/gt.m
@@ -42,26 +42,40 @@ end
 if isequal(A.ndat, B.ndat)
     C = NaN(A.ndat,1);
     for i=1:A.ndat
-        C(i) = compare_vectors(@gt, A.time(i,:), B.time(i,:));
+        C(i) = greaterthan(A.time(i,:), B.time(i,:));
     end
 else
-    if isequal(A.ndat,1) && isequal(B.ndat,1)
-        C = compare_vectors(@gt, A.time, B.time);
-    elseif isequal(A.ndat,1)
+    if isequal(A.ndat,1)
         C = NaN(B.ndat,1);
         for i=1:B.ndat
-            C(i) = compare_vectors(@gt, A.time, B.time(i,:));
+            C(i) = greaterthan(A.time, B.time(i,:));
         end
     elseif isequal(B.ndat,1)
         C = NaN(A.ndat,1);
         for i=1:A.ndat
-            C(i) = compare_vectors(@gt, A.time(i,:), B.time);
+            C(i) = greaterthan(A.time(i,:), B.time);
         end
     else
         C = 0;
     end
 end
 
+
+function c = greaterthan(a,b)
+    if a(1)>b(1)
+        c = 1;
+    else
+        if a(1)<b(1)
+            c = 0;
+        else
+            if a(2)>b(2)
+                c = 1;
+            else
+                c = 0;
+            end
+        end
+    end
+
 %@test:1
 %$ % Define some dates
 %$ date_2 = '1950Q2';
diff --git a/matlab/@dates/le.m b/matlab/@dates/le.m
index 228317228..3056785fe 100644
--- a/matlab/@dates/le.m
+++ b/matlab/@dates/le.m
@@ -43,26 +43,41 @@ if isequal(A.ndat, B.ndat)
     C = (A==B);
     idx = find(C==0);
     for i=1:length(idx)
-        C(idx(i)) = compare_vectors(@lt, A.time(idx(i),:), B.time(idx(i),:));
+        C(idx(i)) = lessorequal(A.time(idx(i),:), B.time(idx(i),:));
     end
 else
-    if isequal(A.ndat,1) && isequal(B.ndat,1)
-        C = compare_vectors(@le, A.time, B.time);
-    elseif isequal(A.ndat,1)
+    if isequal(A.ndat,1)
         C = NaN(B.ndat,1);
         for i=1:B.ndat
-            C(i) = compare_vectors(@le, A.time, B.time(i,:));
+            C(i) = lessorequal(A.time, B.time(i,:));
         end
     elseif isequal(B.ndat,1)
         C = NaN(A.ndat,1);
         for i=1:A.ndat
-            C(i) = compare_vectors(@le, A.time(i,:), B.time);
+            C(i) = lessorequal(A.time(i,:), B.time);
         end
     else
         C = 0;
     end
 end
 
+
+function c = lessorequal(a, b)
+    if a(1)<b(1)
+        c = 1;
+    else
+        if a(1)>b(1)
+            c = 0;
+        else
+            if a(2)<=b(2)
+                c = 1;
+            else
+                c = 0;
+            end
+        end
+    end
+
+
 %@test:1
 %$ % Define some dates
 %$ date_2 = '1950Q2';
@@ -79,7 +94,7 @@ end
 %$ i2 = (d3<=d4);
 %$ i3 = (d4<=d2);
 %$ i4 = (d5<=d4);
-%$ i5 = (d5<=d5); 
+%$ i5 = (d5<=d5);
 %$
 %$ % Check the results.
 %$ t(1) = dyn_assert(i1,1);
diff --git a/matlab/@dates/lt.m b/matlab/@dates/lt.m
index 3707b6861..97bddcff8 100644
--- a/matlab/@dates/lt.m
+++ b/matlab/@dates/lt.m
@@ -42,26 +42,39 @@ end
 if isequal(A.ndat, B.ndat)
     C = NaN(A.ndat,1);
     for i=1:A.ndat
-        C(i) = compare_vectors(@lt, A.time(i,:), B.time(i,:));
+        C(i) = lessthan(A.time(i,:),B.time(i,:));
     end
 else
-    if isequal(A.ndat,1) && isequal(B.ndat,1)
-        C = compare_vectors(@lt, A.time, B.time);
-    elseif isequal(A.ndat,1)
+    if isequal(A.ndat,1)
         C = NaN(B.ndat,1);
         for i=1:B.ndat
-            C(i) = compare_vectors(@lt, A.time, B.time(i,:));
+            C(i) = lessthan(A.time,B.time(i,:));
         end
     elseif isequal(B.ndat,1)
         C = NaN(A.ndat,1);
         for i=1:A.ndat
-            C(i) = compare_vectors(@lt, A.time(i,:), B.time);
+            C(i) =  lessthan(A.time(i,:),B.time);
         end
     else
         C = 0;
     end
 end
 
+function c = lessthan(a,b)
+    if a(1)<b(1)
+        c = 1;
+    else
+        if a(1)>b(1)
+            c = 0;
+        else
+            if a(2)<b(2)
+                c = 1;
+            else
+                c = 0;
+            end
+        end
+    end
+
 %@test:1
 %$ % Define some dates
 %$ date_2 = '1950Q2';
diff --git a/matlab/@dates/subsref.m b/matlab/@dates/subsref.m
index 0cd0cd401..e966fc1e9 100644
--- a/matlab/@dates/subsref.m
+++ b/matlab/@dates/subsref.m
@@ -93,7 +93,7 @@ switch S(1).type
                 elseif isequal(m,1)
                     B.time = [S(1).subs{2}, ones(n,1)];
                 else
-                    error(['dates::subsref: This is a bug!'])
+                    error('dates::subsref: This is a bug!')
                 end
                 B.ndat = rows(B.time);
             elseif isequal(length(S(1).subs),3)
@@ -148,7 +148,7 @@ switch S(1).type
                 elseif isequal(m,1) && isequal(B.freq,1)
                     B.time = [S(1).subs{1}, ones(n,1)];
                 else
-                    error(['dates::subsref: This is a bug!'])
+                    error('dates::subsref: This is a bug!')
                 end
                 B.ndat = rows(B.time);
             else
@@ -267,7 +267,7 @@ end
 %$
 %$ % Define a ranges of dates using qq.
 %$ try
-%$     r1 = qq(1950,1):qq([1950, 3]);
+%$     r1 = qq(1950,1):qq(1950,3);%qq([1950, 3]);
 %$     t(1) = 1;
 %$ catch
 %$     t(1) = 0;
diff --git a/matlab/@dseries/abs.m b/matlab/@dseries/abs.m
index b8d49bce7..72a9520ce 100644
--- a/matlab/@dseries/abs.m
+++ b/matlab/@dseries/abs.m
@@ -40,10 +40,8 @@ function A = abs(B) % --*-- Unitary tests --*--
 
 A = dseries();
 
-A.freq = B.freq;
 A.nobs = B.nobs;
 A.vobs = B.vobs;
-A.init = B.init;
 A.dates = B.dates;
 A.name = cell(A.vobs,1);
 A.tex = cell(A.vobs,1);
diff --git a/matlab/@dseries/align.m b/matlab/@dseries/align.m
index 60c65a947..e785e081d 100644
--- a/matlab/@dseries/align.m
+++ b/matlab/@dseries/align.m
@@ -44,50 +44,50 @@ function [a,b] = align(a, b) % --*-- Unitary tests --*--
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-if ~isequal(a.freq,b.freq)
+if ~isequal(frequency(a),frequency(b))
     error(['dseries::align: ''' inputname(1) ''' and ''' inputname(2) ''' dseries objects must have common frequencies!'])
 end
 
-init = min(a.init,b.init);
+init = min(firstdate(a),firstdate(b));
+last = max(lastdate(a),lastdate(b));
 
-time_range_of_a = a.init:a.init+a.nobs;
-time_range_of_b = b.init:b.init+b.nobs;
-
-last_a = time_range_of_a(a.nobs);
-last_b = time_range_of_b(b.nobs);
-
-common_time_range = intersect(time_range_of_a,time_range_of_b);
-
-if isempty(common_time_range)
+if isempty(intersect(a.dates,b.dates))
     error(['dseries::align: ''' inputname(1) ''' and ''' inputname(2) ''' dseries object must have at least one common date!'])
 end
 
-if a.init<b.init
-    n = b.init-a.init;
+a_init = init;
+b_init = init;
+a_last = last;
+b_last = last;
+
+if firstdate(b)>init
+    n = firstdate(b)-init;
     b.data = [NaN(n,b.vobs); b.data];
     b.nobs = b.nobs+n;
-    b.init = init;
+    b_init = init;
 end
 
-if a.init>b.init
-    n = a.init-b.init;
+if firstdate(a)>init
+    n = firstdate(a)-init;
     a.data = [NaN(n,a.vobs); a.data];
     a.nobs = a.nobs+n;
-    a.init = init;
+    a_init = init;
 end
 
-if last_a>last_b
-    n = last_a-last_b;
+if lastdate(b)<last
+    n = last-lastdate(b);
     b.data = [b.data; NaN(n,b.vobs)];
     b.nobs = b.nobs+n;
-elseif last_a<last_b
-    n = last_b-last_a;
+end
+
+if lastdate(a)<last
+    n = last-lastdate(a);
     a.data = [a.data; NaN(n,a.vobs)];
     a.nobs = a.nobs+n;
 end
 
-a.dates = a.init:a.init+(a.nobs-1);
-b.dates = b.init:b.init+(b.nobs-1);
+a.dates = a_init:a_init+(a.nobs-1);
+b.dates = b_init:b_init+(b.nobs-1);
 
 %@test:1
 %$ % Define a datasets.
@@ -111,7 +111,7 @@ b.dates = b.init:b.init+(b.nobs-1);
 %$ catch
 %$   t(1) = 0;
 %$ end
-%$ 
+%$
 %$ if t(1)
 %$   t(2) = dyn_assert(ts1.nobs,ts2.nobs);
 %$   t(3) = dyn_assert(isequal(ts1.init,ts2.init),1);
diff --git a/matlab/@dseries/baxter_king_filter.m b/matlab/@dseries/baxter_king_filter.m
index 8dce6288a..d019e38c7 100644
--- a/matlab/@dseries/baxter_king_filter.m
+++ b/matlab/@dseries/baxter_king_filter.m
@@ -94,8 +94,8 @@ end
 % Update dseries object.
 ts.data = tmp(K+1:end-K,:);
 ts.nobs = ts.nobs-2*K;
-ts.init = ts.init+K;
-ts.dates = ts.init:ts.init+(ts.nobs-1);
+init = firstdate(ts)+K;
+ts.dates = init:init+(ts.nobs-1);
 
 %@test:1
 %$ plot_flag = 0;
diff --git a/matlab/@dseries/chain.m b/matlab/@dseries/chain.m
index dccf3bb52..2c261758e 100644
--- a/matlab/@dseries/chain.m
+++ b/matlab/@dseries/chain.m
@@ -1,4 +1,4 @@
-function vs = chain(ts,us)
+function vs = chain(ts,us)  % --*-- Unitary tests --*--
 
 % Copyright (C) 2014 Dynare Team
 %
@@ -21,11 +21,11 @@ if ts.vobs-us.vobs
     error(['dseries::chain: dseries objects ' inputname(1) ' and ' inputname(2) ' must have the same number of variables!'])
 end
 
-if ts.freq-us.freq
+if frequency(ts)-frequency(us)
     error(['dseries::chain: dseries objects ' inputname(1) ' and ' inputname(2) ' must have common frequencies!'])
 end
 
-if ts.dates(end)<us.dates(1)
+if lastdate(ts)<firstdate(us)
     error(['dseries::chain: The last date in ' inputname(1) ' (' date2string(ts.dates(end)) ') must not preceed the first date in ' inputname(2) ' (' date2string(us.dates(1)) ')!'])
 end
 
@@ -37,7 +37,7 @@ vs = ts;
 vs.data = [vs.data; bsxfun(@times,CumulatedGrowthFactors,vs.data(end,:))];
 vs.nobs = rows(vs.data);
 
-vs.dates = vs.init:(vs.init+vs.nobs);
+vs.dates = firstdate(vs):firstdate(vs)+vs.nobs;
 
 %@test:1
 %$ try
diff --git a/matlab/@dseries/check.m b/matlab/@dseries/check.m
index 9f07b0119..a190b0dcc 100644
--- a/matlab/@dseries/check.m
+++ b/matlab/@dseries/check.m
@@ -69,7 +69,7 @@ if ~isequal(numel(unique(A.name)),numel(A.name));
     return
 end
 
-if ~isequa(numel(unique(A.tex)),numel(A.tex));
+if ~isequal(numel(unique(A.tex)),numel(A.tex));
     error_flag = 1;
     if nargout>1
         message = ['dseries: The variable tex names in dseries object ''' inputname(1) ''' are not unique!'];
@@ -77,7 +77,7 @@ if ~isequa(numel(unique(A.tex)),numel(A.tex));
     return
 end
 
-if ~isequal(A.dates,A.init:A.init+A.nobs)
+if ~isequal(A.dates,firstdate(A):firstdate(A)+A.nobs)
     error_flag = 1;
     if nargout>1
         message = ['dseries: Wrong definition of the dates member in dseries object ''' inputname(1) '''!'];
diff --git a/matlab/@dseries/cumprod.m b/matlab/@dseries/cumprod.m
index 54a4e33ad..35ed8275e 100644
--- a/matlab/@dseries/cumprod.m
+++ b/matlab/@dseries/cumprod.m
@@ -67,7 +67,7 @@ switch nargin
         B = cumprod(varargin{1});
         t = find(B.dates==varargin{2});
         if isempty(t)
-            if varargin{2}==(B.init-1)
+            if varargin{2}==(firstdate(B)-1)
                 return
             else
                 error(['dseries::cumprod: date ' date2string(varargin{2}) ' is not in the sample!'])
@@ -96,7 +96,7 @@ switch nargin
     B = cumprod(varargin{1});
     t = find(B.dates==varargin{2});
     if isempty(t)
-        if varargin{2}==(B.init-1)
+        if varargin{2}==(firstdate(B)-1)
             B.data = bsxfun(@times,B.data,varargin{3}.data);
             return
         else
diff --git a/matlab/@dseries/cumsum.m b/matlab/@dseries/cumsum.m
index 8dc442041..4c082c4ba 100644
--- a/matlab/@dseries/cumsum.m
+++ b/matlab/@dseries/cumsum.m
@@ -66,7 +66,7 @@ switch nargin
         B = cumsum(varargin{1});
         t = find(B.dates==varargin{2});
         if isempty(t)
-            if varargin{2}==(B.init-1)
+            if varargin{2}==(firstdate(B)-1)
                 return
             else
                 error(['dseries::cumsum: date ' date2string(varargin{2}) ' is not in the sample!'])
@@ -95,7 +95,7 @@ switch nargin
     B = cumsum(varargin{1});
     t = find(B.dates==varargin{2});
     if isempty(t)
-        if varargin{2}==(B.init-1)
+        if varargin{2}==(firstdate(B)-1)
             B.data = bsxfun(@plus,B.data,varargin{3}.data);
             return
         else
diff --git a/matlab/@dseries/dseries.m b/matlab/@dseries/dseries.m
index 0257a8931..433b393d5 100644
--- a/matlab/@dseries/dseries.m
+++ b/matlab/@dseries/dseries.m
@@ -76,18 +76,20 @@ function ts = dseries(varargin) % --*-- Unitary tests --*--
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-ts = struct;
+if nargin>0 && ischar(varargin{1}) && isequal(varargin{1},'initialize')
+    ts = struct;
+    ts.data  = [];
+    ts.nobs  = 0;
+    ts.vobs  = 0;
+    ts.name  = {};
+    ts.tex   = {};
+    ts.dates = dates();
+    ts = class(ts,'dseries');
+    assignin('base','emptydseriesobject',ts);
+    return
+end
 
-ts.data  = [];
-ts.nobs  = 0;
-ts.vobs  = 0;
-ts.name  = {};
-ts.tex   = {};
-ts.freq  = [];
-ts.init  = dates();
-ts.dates = dates();
-
-ts = class(ts,'dseries');
+ts = evalin('base','emptydseriesobject');
 
 switch nargin
   case 0
@@ -100,13 +102,10 @@ switch nargin
             error(['dseries::dseries: Input ' inputname(1) ' (identified as a dates object) must be non empty!'])
           case 1
             % Create an empty dseries object with an initial date.
-            ts.init = varargin{1};
-            ts.freq = varargin{1}.freq;
+            ts.dates = varargin{1};
           otherwise
             % A range of dates is passed to the constructor
             ts.dates = varargin{1};
-            ts.init = varargin{1}(1);
-            ts.freq = varargin{1}.freq;
         end
         return
     elseif ischar(varargin{1})
@@ -136,12 +135,11 @@ switch nargin
         else
             error(['dseries:: I''m not able to load data from ' inputname(1) '!'])
         end
-        ts.init = init;
-        ts.freq = freq;
         ts.data = data;
         ts.name = varlist;
         ts.vobs = length(varlist);
         ts.nobs = size(data,1);
+        ts.dates = init:init+(ts.nobs-1);
         if isempty(tex)
             ts.tex = name2tex(varlist);
         else
@@ -150,10 +148,9 @@ switch nargin
     elseif isnumeric(varargin{1}) && isequal(ndims(varargin{1}),2)
         ts.data = varargin{1};
         [ts.nobs, ts.vobs] = size(ts.data);
-        ts.freq = 1;
-        ts.init = dates(1,1);
         ts.name = default_name(ts.vobs);
         ts.tex = name2tex(ts.name);
+        ts.dates = dates(1,1):dates(1,1)+(ts.nobs-1);
     end
   case {2,3,4}
     a = varargin{1};
@@ -180,20 +177,15 @@ switch nargin
     ts.vobs = size(a,2);
     % Get the first date and set the frequency.
     if isempty(b)
-        ts.freq = 1;
-        ts.init = dates(1,1);
+        init = dates(1,1);
     elseif (isdates(b) && isequal(length(b),1))
-        ts.freq = b.freq;
-        ts.init = b;
+        init = b;
     elseif isdate(b)% Weekly, Monthly, Quaterly or Annual data (string).
-        ts.init = dates(b);
-        ts.freq = ts.init.freq;
+        init = dates(b);
     elseif (isnumeric(b) && isscalar(b) && isint(b)) % Yearly data.
-        ts.freq = 1;
-        ts.init = dates([num2str(b) 'Y']);
+        init = dates([num2str(b) 'Y']);
     elseif isdates(b) % Range of dates
-        ts.freq = b.freq;
-        ts.init = b(1);
+        init = b(1);
         if ts.nobs>1 && ~isequal(b.ndat,ts.nobs)
             message =   'dseries::dseries: If second input is a range, its number of elements must match ';
             message = char(message, '                  the number of rows in the first input, unless the first input');
@@ -218,7 +210,7 @@ switch nargin
     if ~isempty(c)
         if ts.vobs==length(c)
             for i=1:ts.vobs
-                ts.name = vertcat(ts.name, c(i) );
+                ts.name = vertcat(ts.name, c(i));
             end
         else
             error('dseries::dseries: The number of declared names does not match the number of variables!')
@@ -242,7 +234,7 @@ switch nargin
 end
 
 if isempty(ts.dates)
-    ts.dates = ts.init:ts.init+(ts.nobs-1);
+    ts.dates = init:init+(ts.nobs-1);
 end
 
 %@test:1
diff --git a/matlab/@dseries/end.m b/matlab/@dseries/end.m
new file mode 100644
index 000000000..a90869288
--- /dev/null
+++ b/matlab/@dseries/end.m
@@ -0,0 +1,21 @@
+function lastIndex = end(o, k, n)
+
+% Copyright (C) 2014 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+assert(k==1 && n==1, 'dseries::end: Wrong indexing!');
+lastIndex = o.vobs;
\ No newline at end of file
diff --git a/matlab/@dseries/eq.m b/matlab/@dseries/eq.m
index 0f961fb55..7dcd86b2a 100644
--- a/matlab/@dseries/eq.m
+++ b/matlab/@dseries/eq.m
@@ -49,13 +49,13 @@ if ~isequal(A.vobs,B.vobs)
     return
 end
 
-if ~isequal(A.freq,B.freq)
+if ~isequal(frequency(A),frequency(B))
     warning('dseries::eq: Both input arguments should have the same frequencies!')
     C = 0;
     return
 end
 
-if ~isequal(A.init,B.init)
+if ~isequal(firstdate(A),firstdate(B))
     warning('dseries::eq: Both input arguments should have the same initial period!')
     C = 0;
     return
diff --git a/matlab/@dseries/extract.m b/matlab/@dseries/extract.m
index 2c0e0be79..2d61a0b95 100644
--- a/matlab/@dseries/extract.m
+++ b/matlab/@dseries/extract.m
@@ -126,8 +126,6 @@ end
 
 A.data = B.data(:,idVariableName);
 A.dates = B.dates;
-A.init = B.init;
-A.freq = B.freq;
 A.nobs = B.nobs;
 A.vobs = length(idVariableName);
 A.name = B.name(idVariableName);
diff --git a/matlab/@dseries/firstdate.m b/matlab/@dseries/firstdate.m
new file mode 100644
index 000000000..9a7db50ea
--- /dev/null
+++ b/matlab/@dseries/firstdate.m
@@ -0,0 +1,20 @@
+function f = firstdate(o)
+
+% Copyright (C) 2014 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+f = o.dates(1);
\ No newline at end of file
diff --git a/matlab/@dseries/frequency.m b/matlab/@dseries/frequency.m
new file mode 100644
index 000000000..c2befb547
--- /dev/null
+++ b/matlab/@dseries/frequency.m
@@ -0,0 +1,20 @@
+function f = frequency(o)
+
+% Copyright (C) 2014 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+f = o.dates.freq;
\ No newline at end of file
diff --git a/matlab/@dseries/horzcat.m b/matlab/@dseries/horzcat.m
index 40610c3f3..5f05af87f 100644
--- a/matlab/@dseries/horzcat.m
+++ b/matlab/@dseries/horzcat.m
@@ -51,14 +51,21 @@ end
 
 function a = concatenate(b,c)
     [n,message] = common_strings_in_cell_arrays(b.name,c.name);
+    if isempty(b)
+        a = c;
+        return
+    end
+    if isempty(c)
+        a = b;
+        return
+    end
     if n
         error(['dseries::horzcat: I cannot concatenate dseries objects with common variable names (' message ')!'])
     end
-    if ~isequal(b.freq,c.freq)
+    if ~isequal(frequency(b),frequency(c))
         error('dseries::horzcat: All time series objects must have common frequency!')
     else
         a = dseries();
-        a.freq = b.freq;
     end
     d_nobs_flag = 0;
     if ~isequal(b.nobs,c.nobs)
@@ -67,29 +74,26 @@ function a = concatenate(b,c)
         a.nobs = b.nobs;
     end
     d_init_flag = 0;
-    if ~isequal(b.init,c.init)
+    if ~isequal(firstdate(b),firstdate(c))
         d_init_flag = 1;
     end
     a.vobs = b.vobs+c.vobs;
     a.name = vertcat(b.name,c.name);
     a.tex  = vertcat(b.tex,c.tex);
     if ~( d_nobs_flag(1) || d_init_flag(1) )
-        a.init = b.init;
         a.data = [b.data,c.data];
         a.dates = b.dates;
     else
-        if b.init<=c.init
-            a.init = b.init;
-            if b.init<c.init
-                c.data = [NaN(c.init-b.init,c.vobs); c.data];
+        if firstdate(b)<=firstdate(c)
+            if firstdate(b)<firstdate(c)
+                c.data = [NaN(firstdate(c)-firstdate(b),c.vobs); c.data];
             end
         else
-            a.init = c.init;
-            b_first_lines = b.init-c.init;
-            b.data = [NaN(b.init-c.init,b.vobs); b.data];
+            b_first_lines = firstdate(b)-firstdate(c);
+            b.data = [NaN(firstdate(b)-firstdate(c),b.vobs); b.data];
         end
-        b_last_date = b.init+b.nobs;
-        c_last_date = c.init+c.nobs;
+        b_last_date = firstdate(b)+b.nobs;
+        c_last_date = firstdate(c)+c.nobs;
         if b_last_date<c_last_date
             b.data = [b.data; NaN(c_last_date-b_last_date,b.vobs)];
         elseif b_last_date>c_last_date
@@ -289,3 +293,28 @@ function a = concatenate(b,c)
 %$
 %$ T = t;
 %@eof:6
+
+%@test:7
+%$ % Define X
+%$ X = randn(30,2);
+%$
+%$ % Instantiate two time series objects.
+%$ ts1 = dseries();
+%$ ts2 = dseries(randn(30,2),'1950Q2');
+%$
+%$ % Call the tested method.
+%$ try
+%$   ts3 = [ts1,ts2];
+%$   t = 1;
+%$ catch
+%$   t = 0;
+%$ end
+%$
+%$ if t(1)
+%$   t(2) = dyn_assert(ts3.freq,4);
+%$   t(3) = dyn_assert(ts3.data,X);
+%$   t(4) = dyn_assert(isequal(ts3.dates(1),dates('1950Q2')),1);
+%$ end
+%$
+%$ T = t;
+%@eof:7
diff --git a/matlab/@dseries/insert.m b/matlab/@dseries/insert.m
index 4b8aa06b9..128c72dff 100644
--- a/matlab/@dseries/insert.m
+++ b/matlab/@dseries/insert.m
@@ -52,7 +52,7 @@ if n
     error(['dseries::insert: Variable(s) ' message ' already exist in ''' inputname(1) '''!'])
 end
 
-if ~isequal(ts.freq,us.freq)
+if ~isequal(frequency(ts),frequency(us))
     error(['dseries::insert: ''' inputname(1) ''' and ''' inputname(2) ''' dseries objects must have common frequencies!'])
 end
 
diff --git a/matlab/@dseries/isequal.m b/matlab/@dseries/isequal.m
index 50b93b097..003f6903c 100644
--- a/matlab/@dseries/isequal.m
+++ b/matlab/@dseries/isequal.m
@@ -44,12 +44,12 @@ if ~isequal(A.vobs,B.vobs)
     return
 end
 
-if ~isequal(A.freq,B.freq)
+if ~isequal(frequency(A),frequency(B))
     C = 0;
     return
 end
 
-if ~isequal(A.init,B.init)
+if ~isequal(A.dates,B.dates)
     C = 0;
     return
 end
diff --git a/matlab/@dseries/lastdate.m b/matlab/@dseries/lastdate.m
new file mode 100644
index 000000000..972a70e18
--- /dev/null
+++ b/matlab/@dseries/lastdate.m
@@ -0,0 +1,20 @@
+function l = lastdate(o)
+
+% Copyright (C) 2014 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+l = o.dates(end);
\ No newline at end of file
diff --git a/matlab/@dseries/merge.m b/matlab/@dseries/merge.m
index 52fe8ac48..0f02ceea8 100644
--- a/matlab/@dseries/merge.m
+++ b/matlab/@dseries/merge.m
@@ -44,12 +44,11 @@ if ~isdseries(C)
     error('dseries::merge: Both inputs must be dseries objects!')
 end
 
-if ~isequal(B.freq,C.freq)
+if ~isequal(frequency(B),frequency(C))
     error(['dseries::merge: Cannot merge ' inputname(1) ' and ' inputname(2) ' (frequencies are different)!'])
 end
 
 A = dseries();
-A.freq = B.freq;
 [A.name, IBC, junk] = unique([B.name; C.name], 'last');
 tex = [B.tex; C.tex];
 A.tex = tex(IBC); 
@@ -59,8 +58,8 @@ if B.nobs == 0
     A = C;
 elseif C.nobs == 0
     A = B;
-elseif B.init >= C.init
-    diff = B.init - C.init;
+elseif firstdate(B) >= firstdate(C)
+    diff = firstdate(B) - firstdate(C);
     A.nobs = max(B.nobs + diff, C.nobs);
     A.data = NaN(A.nobs, A.vobs);
     Z1 = [NaN(diff,B.vobs);B.data];
@@ -68,14 +67,14 @@ elseif B.init >= C.init
         Z1 = [Z1; NaN(A.nobs-(B.nobs + diff),B.vobs)];
     end;
     Z2 = C.data;
-    if A.nobs > C.nobs 
+    if A.nobs > C.nobs
         Z2 = [Z2; NaN(A.nobs - C.nobs,C.vobs)];
     end;
     Z = [Z1 Z2];
     A.data = Z(:,IBC);
-    A.init = C.init;
+    A_init = firstdate(C);
 else
-    diff = C.init - B.init;
+    diff = firstdate(C) - firstdate(B);
     A.nobs = max(C.nobs + diff, B.nobs);
     A.data = NaN(A.nobs, A.vobs);
     Z1 = [NaN(diff,C.vobs);C.data];
@@ -88,10 +87,10 @@ else
     end;
     Z = [Z2 Z1];
     A.data = Z(:,IBC);
-    A.init = B.init;
+    A_init = B.init;
 end
 
-A.dates = A.init:A.init+(A.nobs-1);
+A.dates = A_init:A_init+(A.nobs-1);
 
 %@test:1
 %$ % Define a datasets.
diff --git a/matlab/@dseries/minus.m b/matlab/@dseries/minus.m
index 6f7c47dbc..f563485bc 100644
--- a/matlab/@dseries/minus.m
+++ b/matlab/@dseries/minus.m
@@ -73,11 +73,11 @@ else
     end
 end
 
-if ~isequal(B.freq,C.freq)
+if ~isequal(frequency(B),frequency(C))
     error(['dseries::plus: Cannot substract ' inputname(1) ' and ' inputname(2) ' (frequencies are different)!'])
 end
 
-if ~isequal(B.nobs,C.nobs) || ~isequal(B.init,C.init)
+if ~isequal(B.nobs,C.nobs) || ~isequal(firstdate(B),firstdate(C))
     [B, C] = align(B, C);
 end
 
@@ -93,8 +93,6 @@ end
 
 A = dseries();
 
-A.freq = B.freq;
-A.init = B.init;
 A.dates = B.dates;
 A.nobs = max(B.nobs,C.nobs);
 A.vobs = max(B.vobs,C.vobs);
diff --git a/matlab/@dseries/mpower.m b/matlab/@dseries/mpower.m
index 89b0b1049..06cc5af73 100644
--- a/matlab/@dseries/mpower.m
+++ b/matlab/@dseries/mpower.m
@@ -40,10 +40,26 @@ function A = mpower(B,C) % --*-- Unitary tests --*--
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-if isdseries(B) && isnumeric(C) && isreal(C) && isscalar(C)
+if isnumeric(B) && isvector(B) && length(B)>1
+    if ~isdseries(C)
+        error('dseries::mpower: Second input argument must be a dseries object!')
+    end
+    A = C;
+    A.data = bsxfun(@power,C.data,B);
+    return;
+end
+
+if isnumeric(C) && isvector(C) && length(C)>1
+    if ~isdseries(B)
+        error('dseries::mpower: First input argument must be a dseries object!')
+    end
+    A = B;
+    A.data = bsxfun(@power,B.data,C);
+    return
+end
+
+if isdseries(B) && isnumeric(C) && isreal(C) &&  isscalar(C)
     A = dseries();
-    A.freq = B.freq;
-    A.init = B.init;
     A.dates = B.dates;
     A.nobs = B.nobs;
     A.vobs = B.vobs;
@@ -58,10 +74,8 @@ if isdseries(B) && isnumeric(C) && isreal(C) && isscalar(C)
 end
 
 if isdseries(B) && isdseries(C)
-    if isequal(B.nobs,C.nobs) && isequal(B.vobs,C.vobs) && isequal(B.freq,C.freq)
+    if isequal(B.nobs,C.nobs) && isequal(B.vobs,C.vobs) && isequal(frequency(B),frequency(C))
         A = dseries();
-        A.freq = B.freq;
-        A.init = B.init;
         A.dates = B.dates;
         A.nobs = B.nobs;
         A.vobs = B.vobs;
@@ -133,4 +147,26 @@ error(['dseries::mpower: Wrong calling sequence!'])
 %$    t(6) = dyn_assert(ts3.tex,{'A1^2';'A2^2'});
 %$ end
 %$ T = all(t);
-%@eof:2
\ No newline at end of file
+%@eof:2
+
+%@test:3
+%$ % Define a dseries object
+%$ ts1=dseries([1 1;2 2;3 3], '1999y', {'MyVar1','MyVar2'});
+%$
+%$ % Use the power
+%$ try
+%$    ts2 = ts1^transpose(1:3);
+%$    t = 1;
+%$ catch
+%$    t = 0;
+%$ end
+%$
+%$ if t(1)
+%$    t(2) = dyn_assert(ts2.vobs,2);
+%$    t(3) = dyn_assert(ts2.nobs,3);
+%$    t(4) = dyn_assert(ts2.data,bsxfun(@power,ts1.data,transpose(1:3)),1e-15);
+%$    t(5) = dyn_assert(ts2.name,{'MyVar1';'MyVar2'});
+%$    t(6) = dyn_assert(ts2.tex,{'MyVar1';'MyVar2'});
+%$ end
+%$ T = all(t);
+%@eof:3
\ No newline at end of file
diff --git a/matlab/@dseries/mrdivide.m b/matlab/@dseries/mrdivide.m
index e42b47aa9..47d4c2127 100644
--- a/matlab/@dseries/mrdivide.m
+++ b/matlab/@dseries/mrdivide.m
@@ -74,15 +74,13 @@ if isdseries(B) && isdseries(C)
             idC = 1:C.vobs;
         end
     end
-    if ~isequal(B.freq,C.freq)
+    if ~isequal(frequency(B),frequency(C))
         error(['dseries::times: Cannot divide ' inputname(1) ' and ' inputname(2) ' (frequencies are different)!'])
     end
-    if ~isequal(B.nobs,C.nobs) || ~isequal(B.init,C.init)
+    if ~isequal(B.nobs,C.nobs) || ~isequal(firstdate(B),firstdate(C))
         [B, C] = align(B, C);
     end
     A = dseries();
-    A.freq = B.freq;
-    A.init = B.init;
     A.dates = B.dates;
     A.nobs = max(B.nobs,C.nobs);
     A.vobs = max(B.vobs,C.vobs);
diff --git a/matlab/@dseries/mtimes.m b/matlab/@dseries/mtimes.m
index 79f2019c2..719687efd 100644
--- a/matlab/@dseries/mtimes.m
+++ b/matlab/@dseries/mtimes.m
@@ -74,15 +74,13 @@ if isdseries(B) && isdseries(C)
             idC = 1:C.vobs;
         end
     end
-    if ~isequal(B.freq,C.freq)
+    if ~isequal(frequency(B),frequency(C))
         error(['dseries::times: Cannot multiply ' inputname(1) ' and ' inputname(2) ' (frequencies are different)!'])
     end
-    if ~isequal(B.nobs,C.nobs) || ~isequal(B.init,C.init)
+    if ~isequal(B.nobs,C.nobs) || ~isequal(firstdate(B),firstdate(C))
         [B, C] = align(B, C);
     end
     A = dseries();
-    A.freq = B.freq;
-    A.init = B.init;
     A.dates = B.dates;
     A.nobs = max(B.nobs,C.nobs);
     A.vobs = max(B.vobs,C.vobs);
diff --git a/matlab/@dseries/ne.m b/matlab/@dseries/ne.m
index 53c98867c..b8056d82f 100644
--- a/matlab/@dseries/ne.m
+++ b/matlab/@dseries/ne.m
@@ -38,35 +38,35 @@ if ~(isdseries(A) && isdseries(B))
 end
 
 if ~isequal(A.nobs,B.nobs)
-    warning('dseries::eq: Both input arguments should have the same number of observations!')
+    warning('dseries::ne: Both input arguments should have the same number of observations!')
     C = 1;
     return
 end
 
 if ~isequal(A.vobs,B.vobs)
-    warning('dseries::eq: Both input arguments should have the same number of observations!')
+    warning('dseries::ne: Both input arguments should have the same number of observations!')
     C = 1;
     return
 end
 
-if ~isequal(A.freq,B.freq)
-    warning('dseries::eq: Both input arguments should have the same frequencies!')
+if ~isequal(frequency(A),frequency(B))
+    warning('dseries::ne: Both input arguments should have the same frequencies!')
     C = 1;
     return
 end
 
-if ~isequal(A.init,B.init)
-    warning('dseries::eq: Both input arguments should have the same initial period!')
+if ~isequal(firstdate(A),firstdate(B))
+    warning('dseries::ne: Both input arguments should have the same initial period!')
     C = 1;
     return
 end
 
 if ~isequal(A.name,B.name)
-    warning('dseries::eq: Both input arguments do not have the same variables!')
+    warning('dseries::ne: Both input arguments do not have the same variables!')
 end
 
 if ~isequal(A.tex,B.tex)
-    warning('dseries::eq: Both input arguments do not have the same tex names!')
+    warning('dseries::ne: Both input arguments do not have the same tex names!')
 end
 
 C = ne(A.data, B.data);
diff --git a/matlab/@dseries/plus.m b/matlab/@dseries/plus.m
index b9d20a798..be116601c 100644
--- a/matlab/@dseries/plus.m
+++ b/matlab/@dseries/plus.m
@@ -73,11 +73,11 @@ else
     end
 end
 
-if ~isequal(B.freq,C.freq)
+if ~isequal(frequency(B),frequency(C))
     error(['dseries::plus: Cannot add ' inputname(1) ' and ' inputname(2) ' (frequencies are different)!'])
 end
 
-if ~isequal(B.nobs,C.nobs) || ~isequal(B.init,C.init)
+if ~isequal(B.nobs,C.nobs) || ~isequal(firstdate(B),firstdate(C))
     [B, C] = align(B, C);
 end
 
@@ -93,8 +93,9 @@ end
 
 A = dseries();
 
-A.freq = B.freq;
-A.init = B.init;
+A.dates = B.dates;
+%A.freq = B.freq;
+%A.init = B.init;
 A.nobs = max(B.nobs,C.nobs);
 A.vobs = max(B.vobs,C.vobs);
 A.name = cell(A.vobs,1);
@@ -104,7 +105,7 @@ for i=1:A.vobs
     A.tex(i) = {['(' B.tex{idB(i)} '+' C.tex{idC(i)} ')']};
 end
 A.data = bsxfun(@plus,B.data,C.data);
-A.dates = A.init:A.init+(A.nobs-1);
+%A.dates = A.init:A.init+(A.nobs-1);
 
 %@test:1
 %$ % Define a datasets.
diff --git a/matlab/@dseries/qdiff.m b/matlab/@dseries/qdiff.m
index b1683f1d6..f57410e31 100644
--- a/matlab/@dseries/qdiff.m
+++ b/matlab/@dseries/qdiff.m
@@ -41,7 +41,7 @@ function us = qdiff(ts) % --*-- Unitary tests --*--
 
 us = ts;
 
-switch ts.freq
+switch frequency(ts)
   case 1
     error('dseries::qgrowth: I cannot compute quaterly differences from yearly data!')
   case 4
diff --git a/matlab/@dseries/qgrowth.m b/matlab/@dseries/qgrowth.m
index f86233abe..796a889bf 100644
--- a/matlab/@dseries/qgrowth.m
+++ b/matlab/@dseries/qgrowth.m
@@ -41,7 +41,7 @@ function us = qgrowth(ts) % --*-- Unitary tests --*--
 
 us = ts;
 
-switch ts.freq
+switch frequency(ts)
   case 1
     error('dseries::qgrowth: I cannot compute quaterly growth rates from yearly data!')
   case 4
diff --git a/matlab/@dseries/remove.m b/matlab/@dseries/remove.m
new file mode 100644
index 000000000..db9690cfd
--- /dev/null
+++ b/matlab/@dseries/remove.m
@@ -0,0 +1,98 @@
+function ts = remove(ts,a) % --*-- Unitary tests --*--
+
+% Removes a variable from a dseries object (alias for the pop method).
+
+%@info:
+%! @deftypefn {Function File} {@var{ts} =} pop (@var{ts}, @var{a})
+%! @anchor{dseries/pop}
+%! @sp 1
+%! Remove method for the dseries class. Removes a variable from a dseries object.
+%! @sp 2
+%! @strong{Inputs}
+%! @sp 1
+%! @table @ @var
+%! @item ts
+%! Object instantiated by @ref{dseries}.
+%! @item a
+%! String, name of the variable to be removed.
+%! @end table
+%! @sp 2
+%! @strong{Outputs}
+%! @sp 1
+%! @table @ @var
+%! @item ts
+%! Object instantiated by @ref{dseries}, without variable (@var{a}).
+%! @end table
+%! @end deftypefn
+%@eod:
+
+% Copyright (C) 2014 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+ts = pop(ts, a);
+
+%@test:1
+%$ % Define a datasets.
+%$ A = rand(10,3);
+%$
+%$ % Define names
+%$ A_name = {'A1';'A2';'A3'};
+%$
+%$ t = zeros(4,1);
+%$
+%$ % Instantiate a time series object.
+%$ try
+%$    ts1 = dseries(A,[],A_name,[]);
+%$    ts2 = remove(ts1,'A2');
+%$    t(1) = 1;
+%$ catch
+%$    t = 0;
+%$ end
+%$
+%$ if length(t)>1
+%$    t(2) = dyn_assert(ts2.vobs,2);
+%$    t(3) = dyn_assert(ts2.nobs,10);
+%$    t(4) = dyn_assert(ts2.data,[A(:,1), A(:,3)],1e-15);
+%$ end
+%$ T = all(t);
+%@eof:1
+
+%@test:2
+%$ % Define a datasets.
+%$ A = rand(10,3);
+%$
+%$ % Define names
+%$ A_name = {'A1';'A2';'A3'};
+%$
+%$ t = zeros(4,1);
+%$
+%$ % Instantiate a time series object.
+%$ try
+%$    ts1 = dseries(A,[],A_name,[]);
+%$    ts2 = ts1.remove('A2');
+%$    t(1) = 1;
+%$ catch
+%$    t = 0;
+%$ end
+%$
+%$ if length(t)>1
+%$    t(2) = dyn_assert(ts2.vobs,2);
+%$    t(3) = dyn_assert(ts2.nobs,10);
+%$    t(4) = dyn_assert(ts2.data,[A(:,1), A(:,3)],1e-15);
+%$ end
+%$ T = all(t);
+%@eof:2
\ No newline at end of file
diff --git a/matlab/@dseries/save.m b/matlab/@dseries/save.m
index cad2b2e8a..2f3ff7c49 100644
--- a/matlab/@dseries/save.m
+++ b/matlab/@dseries/save.m
@@ -35,8 +35,8 @@ switch format
     fid = fopen([basename, '.m'],'w');
     fprintf(fid,'%% File created on %s.\n',datestr(now));
     fprintf(fid,'\n');
-    fprintf(fid,'FREQ__ = %s;\n',num2str(A.freq));
-    fprintf(fid,'INIT__ = '' %s'';\n',date2string(A.init));
+    fprintf(fid,'FREQ__ = %s;\n',num2str(frequency(A)));
+    fprintf(fid,'INIT__ = '' %s'';\n',date2string(firstdate(A)));
     fprintf(fid,'\n');
     fprintf(fid,'NAMES__ = {');
     for i=1:A.vobs
@@ -61,8 +61,8 @@ switch format
     end
     fclose(fid);
   case 'mat'
-    FREQ__ = A.freq;
-    INIT__ = date2string(A.init);
+    FREQ__ = frequency(A);
+    INIT__ = date2string(firstdate(A));
     NAMES__ = A.name;
     TEX__ = A.tex;
     str = [];
diff --git a/matlab/@dseries/subsasgn.m b/matlab/@dseries/subsasgn.m
index 17fc472cd..3937fd092 100644
--- a/matlab/@dseries/subsasgn.m
+++ b/matlab/@dseries/subsasgn.m
@@ -25,7 +25,7 @@ function A = subsasgn(A,S,B) % --*-- Unitary tests --*--
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-merge_dseries_objects = 1;    
+merge_dseries_objects = 1;
 
 switch length(S)
     case 1
@@ -77,6 +77,12 @@ switch length(S)
                   end
               end
           end
+          if isempty(B)
+              for i=1:length(S(1).subs)
+                  A = remove(A,S(1).subs{i});
+              end
+              return
+          end
           if ~isequal(length(S(1).subs),B.vobs)
               error('dseries::subsasgn: Wrong syntax!')
           end
@@ -99,20 +105,14 @@ switch length(S)
           end
         case '.'
           if isequal(S(1).subs,'init') && isdates(B) && isequal(length(B),1)
-              % Overwrite the init member...
-              A.init = B;
-              % ... and update freq and time members.
-              A.freq = A.init.freq;
-              A.dates = A.init:A.init+(A.nobs-1);
+              % Change the initial date (update dates member)
+              A.dates = B:B+(A.nobs-1);
               return
           elseif isequal(S(1).subs,'dates') && isdates(B)
-              % Overwrite the time member...
+              % Overwrite the dates member
               A.dates = B;
-              % ... and update the freq and init members.
-              A.init = B(1);
-              A.freq = A.init.freq;
               return
-          elseif ismember(S(1).subs,{'freq','nobs','vobs','data','name','tex'})
+          elseif ismember(S(1).subs,{'nobs','vobs','data','name','tex'})
               error(['dseries::subsasgn: You cannot overwrite ' S(1).subs ' member!'])
           elseif ~isequal(S(1).subs,B.name)
               % Single variable selection.
@@ -211,7 +211,7 @@ switch length(S)
         else
             sA = extract(A,S(1).subs);
         end
-        if (isdseries(B) && isequal(sA.vobs,B.vobs)) || (isnumeric(B) && isequal(sA.vobs,columns(B))) || (isnumeric(B) && isequal(columns(B),1)) 
+        if (isdseries(B) && isequal(sA.vobs,B.vobs)) || (isnumeric(B) && isequal(sA.vobs,columns(B))) || (isnumeric(B) && isequal(columns(B),1))
             if isdates(S(2).subs{1})
                 [junk, tdx] = intersect(sA.dates.time,S(2).subs{1}.time,'rows');
                 if isdseries(B)
@@ -278,8 +278,8 @@ end
 %$     t(1) = 1;
 %$ catch
 %$     t(1) = 0;
-%$ end 
-%$ 
+%$ end
+%$
 %$ % Instantiate a time series object.
 %$ if t(1)
 %$    t(2) = dyn_assert(ts1.vobs,3);
@@ -440,8 +440,8 @@ end
 %$     t(1) = 1;
 %$ catch
 %$     t(1) = 0;
-%$ end 
-%$ 
+%$ end
+%$
 %$ % Instantiate a time series object.
 %$ if t(1)
 %$    t(2) = dyn_assert(ts1.vobs,4);
@@ -469,8 +469,8 @@ end
 %$     t(1) = 1;
 %$ catch
 %$     t(1) = 0;
-%$ end 
-%$ 
+%$ end
+%$
 %$ % Instantiate a time series object.
 %$ if t(1)
 %$    t(2) = dyn_assert(ts1.vobs,4);
@@ -499,8 +499,8 @@ end
 %$     t(1) = 1;
 %$ catch
 %$     t(1) = 0;
-%$ end 
-%$ 
+%$ end
+%$
 %$ % Instantiate a time series object.
 %$ if t(1)
 %$    t(2) = dyn_assert(ts1.vobs,4);
@@ -845,7 +845,7 @@ end
 %$ T = all(t);
 %@eof:21
 
-%@test:21
+%@test:22
 %$ % Define a datasets.
 %$ A = rand(1,3);
 %$
@@ -867,4 +867,28 @@ end
 %$    t(4) = dyn_assert(ts.data,repmat(A,4,1),1e-15);
 %$ end
 %$ T = all(t);
-%@eof:21
\ No newline at end of file
+%@eof:22
+
+%@test:23
+%$ % Instantiate a dseries object.
+%$ ts0 = dseries(randn(10,6), '1999y');
+%$
+%$ % Try to remove Variable_2 and Variable_3
+%$ try
+%$     ts0{'Variable_@2,3@'} = [];
+%$     t(1) = 1;
+%$ catch
+%$     t(1) = 0;
+%$ end
+%$
+%$ if t(1)
+%$    % Try to display Variable_2 and Variable_3 again (should fail because already removed)
+%$    try
+%$       ts0{'Variable_@2,3@'};
+%$       t(2) = 0;
+%$    catch
+%$       t(2) = 1;
+%$    end
+%$ end
+%$ T = all(t);
+%@eof:23
\ No newline at end of file
diff --git a/matlab/@dseries/subsref.m b/matlab/@dseries/subsref.m
index ab3b08dc5..d2ecdaa90 100644
--- a/matlab/@dseries/subsref.m
+++ b/matlab/@dseries/subsref.m
@@ -65,7 +65,7 @@ function B = subsref(A, S) % --*-- Unitary tests --*--
 switch S(1).type
   case '.'
     switch S(1).subs
-      case {'data','nobs','vobs','name','tex','freq','dates','init'}        % Public members.
+      case {'data','nobs','vobs','name','tex','dates'}        % Public members.
         if length(S)>1 && isequal(S(2).type,'()') && isempty(S(2).subs)
             error(['dseries::subsref: ' S(1).subs ' is not a method but a member!'])
         end
@@ -75,6 +75,15 @@ switch S(1).type
         if length(S)>1 && isequal(S(2).type,'()') && isempty(S(2).subs)
             S = shiftS(S,1);
         end
+      case 'init'
+        % Returns a dates object (first date).
+        B = A.dates(1);
+      case 'last'
+        % Returns a dates object (last date).
+        B = A.dates(end);
+      case 'freq'
+        % Returns an integer characterizing the data frequency (1, 4, 12 or 52)
+        B = A.dates.freq;
       case {'lag','lead','hptrend','hpcycle','chain'} % Methods with less than two arguments.
         if length(S)>1 && isequal(S(2).type,'()')
             if isempty(S(2).subs)
@@ -90,7 +99,7 @@ switch S(1).type
         else
             B = feval(S(1).subs,A);
         end
-      case {'cumsum','insert','pop','cumprod'} % Methods with less than three argument.
+      case {'cumsum','insert','pop','cumprod','remove'} % Methods with less than three argument.
         if length(S)>1 && isequal(S(2).type,'()')
             if isempty(S(2).subs)
                 B = feval(S(1).subs,A);
@@ -170,8 +179,6 @@ switch S(1).type
             B.tex  = deblank(A.tex(ndx,:));
             B.nobs = A.nobs;
             B.vobs = 1;
-            B.freq = A.freq;
-            B.init = A.init;
             B.dates = A.dates;
         else
             error('dseries::subsref: Unknown public method, public member or variable!')
@@ -227,29 +234,9 @@ switch S(1).type
         B.tex  = A.tex;
         B.nobs = length(tdx);
         B.vobs = A.vobs;
-        B.freq = A.freq;
-        B.init = A.init+(tdx(1)-1);
         B.dates = A.dates(tdx);
     elseif isvector(S(1).subs{1}) && all(isint(S(1).subs{1}))
-        % Extract a subsample using a vector of integers (observation index).
-        % Note that this does not work if S(1).subs is an integer scalar... In which case S(1).subs is interpreted as a lead/lag operator (as in the Dynare syntax).
-        % To extract one observation, a dates with one element input must be used.
-        if all(S(1).subs{1}>0) && all(S(1).subs{1}<=A.nobs)
-            if size(A.data,2)>1
-                S(1).subs = [S(1).subs, ':'];
-            end
-            B = dseries();
-            B.data = builtin('subsref', A.data, S(1));
-            B.nobs = size(B.data,1);
-            B.vobs = A.vobs;
-            B.freq = A.freq;
-            B.dates = A.dates(S(1).subs{1});
-            B.init = B.dates(1);
-            B.name = A.name;
-            B.tex  = A.tex;
-        else
-            error('dseries::subsref: Indices are out of bounds!')
-        end
+        error('dseries::subsref: It is not possible to select observations with a vector of integers. You have to index with a dates object instead!');
     else
         error('dseries::subsref: I have no idea of what you are trying to do!')
     end
@@ -267,8 +254,6 @@ switch S(1).type
         B.tex  = A.tex(idx);
         B.nobs = A.nobs;
         B.vobs = length(idx);
-        B.freq = A.freq;
-        B.init = A.init;
         B.dates = A.dates;
     else
         error('dseries::subsref: What the Hell are you tryin'' to do?!')
@@ -293,7 +278,7 @@ end
 %$ ts1 = dseries(A,[],A_name,[]);
 %$
 %$ % Call the tested method.
-%$ a = ts1(2:9);
+%$ a = ts1(ts1.dates(2:9));
 %$
 %$ % Expected results.
 %$ e.data = [transpose(2:9),2*transpose(2:9)];
@@ -653,7 +638,7 @@ end
 %@test:15
 %$ try
 %$     ds = dseries(transpose(1:5));
-%$     ts = ds(2:3);
+%$     ts = ds(ds.dates(2:3));
 %$     t(1) = 1;
 %$ catch
 %$     t(1) = 0;
@@ -670,7 +655,7 @@ end
 %@test:16
 %$ try
 %$     ds = dseries(transpose(1:5));
-%$     ts = ds(2:6);
+%$     ts = ds(ds.dates(2:6));
 %$     t(1) = 0;
 %$ catch
 %$     t(1) = 1;
diff --git a/matlab/@dseries/uminus.m b/matlab/@dseries/uminus.m
index bcccca70a..e8214cd59 100644
--- a/matlab/@dseries/uminus.m
+++ b/matlab/@dseries/uminus.m
@@ -42,10 +42,8 @@ function A = uminus(B) % --*-- Unitary tests --*--
 
 A = dseries();
 
-A.freq = B.freq;
 A.nobs = B.nobs;
 A.vobs = B.vobs;
-A.init = B.init;
 A.dates = B.dates;
 A.name = cell(A.vobs,1);
 A.tex = cell(A.vobs,1);
diff --git a/matlab/@dseries/vertcat.m b/matlab/@dseries/vertcat.m
index 9b286945d..4b4a680fe 100644
--- a/matlab/@dseries/vertcat.m
+++ b/matlab/@dseries/vertcat.m
@@ -60,7 +60,7 @@ end
 
 function d = vertcat_(b, c)
     d = NaN;
-    if ~isequal(b.freq, c.freq)
+    if ~isequal(frequency(b), frequency(c))
         error('dseries::vertcat: Frequencies must be common!')
     end
     if ~isequal(b.vobs, c.vobs)
diff --git a/matlab/@dseries/ydiff.m b/matlab/@dseries/ydiff.m
index 4e20608e6..49b94dc6f 100644
--- a/matlab/@dseries/ydiff.m
+++ b/matlab/@dseries/ydiff.m
@@ -41,7 +41,7 @@ function us = ydiff(ts) % --*-- Unitary tests --*--
 
 us = ts;
 
-switch ts.freq
+switch frequency(ts)
   case 1
     us.data(2:end,:) = ts.data(2:end,:)-ts.data(1:end-1,:);
     us.data(1,:) = NaN;
@@ -55,21 +55,21 @@ switch ts.freq
     for i = 1:ts.vobs
         us.name(i) = {['ydiff(' us.name{i} ')']};
         us.tex(i) = {['\Delta_4 ' us.tex{i}]};
-    end  
+    end
   case 12
     us.data(13:end,:) = ts.data(13:end,:)-ts.data(1:end-12,:);
     us.data(1:12,:) = NaN;
     for i = 1:ts.vobs
         us.name(i) = {['ydiff(' us.name{i} ')']};
         us.tex(i) = {['\Delta_{12} ' us.tex{i}]};
-    end    
+    end
   case 52
     us.data(53:end,:) = ts.data(53:end,:)-ts.data(1:end-52,:);
     us.data(1:52,:) = NaN;
     for i = 1:ts.vobs
         us.name(i) = {['ydiff(' us.name{i} ')']};
         us.tex(i) = {['\Delta_{52} ' us.tex{i}]};
-    end    
+    end
   otherwise
     error(['dseries::ygrowth: object ' inputname(1) ' has unknown frequency']);
 end
diff --git a/matlab/@dseries/ygrowth.m b/matlab/@dseries/ygrowth.m
index 8b26d5252..ba988ec2e 100644
--- a/matlab/@dseries/ygrowth.m
+++ b/matlab/@dseries/ygrowth.m
@@ -41,7 +41,7 @@ function us = ygrowth(ts) % --*-- Unitary tests --*--
 
 us = ts;
 
-switch ts.freq
+switch frequency(ts)
   case 1
     us.data(2:end,:) = ts.data(2:end,:)./ts.data(1:end-1,:) - 1;
     us.data(1,:) = NaN;
@@ -55,7 +55,7 @@ switch ts.freq
     for i = 1:ts.vobs
         us.name(i) = {['ygrowth(' us.name{i} ')']};
         us.tex(i) = {['\delta_4 ' us.tex{i}]};
-    end  
+    end
   case 12
     us.data(13:end,:) = ts.data(13:end,:)./ts.data(1:end-12,:) - 1;
     us.data(1:12,:) = NaN;
diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m
index 9701f3a88..0202f7bfc 100644
--- a/matlab/dynare_config.m
+++ b/matlab/dynare_config.m
@@ -63,6 +63,7 @@ addpath([dynareroot '/utilities/tests/'])
 addpath([dynareroot '/utilities/dates/'])
 addpath([dynareroot '/utilities/dataset/'])
 addpath([dynareroot '/utilities/general/'])
+addpath([dynareroot '/utilities/dseries/'])
 addpath([dynareroot '/reports/'])
 
 % For functions that exist only under some Octave versions
@@ -249,4 +250,8 @@ if verbose
     skipline()
 end
 
+% Save empty dates and dseries objects (necessary if a mod file is not preprocessed).
+dates('initialize');
+dseries('initialize');
+
 cd(origin);
diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index ec2f2c946..f5f166d41 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -64,7 +64,7 @@ if ~isempty(i)
     disp('   i) if all parameters occurring in these equations are defined')
     disp('  ii) that no division by an endogenous variable initialized to 0 occurs') 
     info = 1;
-    x = NaN;
+    x = NaN(size(fvec));
     return;
 end
 
diff --git a/matlab/reports/@graph/graph.m b/matlab/reports/@graph/graph.m
index 687304397..2d9f97560 100644
--- a/matlab/reports/@graph/graph.m
+++ b/matlab/reports/@graph/graph.m
@@ -69,6 +69,10 @@ o.xTickLabels = {};
 o.xTickLabelRotation = 0;
 o.xTickLabelAnchor = 'east';
 
+o.yTickLabelPrecision = 0;
+o.yTickLabelFixed = true;
+o.yTickLabelZeroFill = true;
+
 o.width = 6;
 o.height = 4.5;
 
@@ -125,6 +129,9 @@ assert(isfloat(o.width), '@graph.graph: o.width must be a real number');
 assert(isfloat(o.height), '@graph.graph: o.height must be a real number');
 assert(isfloat(o.xTickLabelRotation), '@graph.graph: o.xTickLabelRotation must be a real number');
 assert(ischar(o.xTickLabelAnchor), '@graph.graph: xTickLabelAnchor must be a string');
+assert(isint(o.yTickLabelPrecision), '@graph.graph: o.yTickLabelPrecision must be an integer');
+assert(islogical(o.yTickLabelFixed), '@graph.graph: yTickLabelFixed must be either true or false');
+assert(islogical(o.yTickLabelZeroFill), '@graph.graph: yTickLabelZeroFill must be either true or false');
 
 valid_shadeColor = {'red', 'green', 'blue', 'cyan ', 'magenta', 'yellow', ...
                     'black', 'gray', 'darkgray', 'lightgray', 'brown', ...
diff --git a/matlab/reports/@graph/writeGraphFile.m b/matlab/reports/@graph/writeGraphFile.m
index 2b706a8fe..bfa480245 100644
--- a/matlab/reports/@graph/writeGraphFile.m
+++ b/matlab/reports/@graph/writeGraphFile.m
@@ -103,7 +103,15 @@ for i = 1:xlen
         fprintf(fid,',');
     end
 end
-fprintf(fid, '},\nx tick label style={rotate=%f', o.xTickLabelRotation);
+fprintf(fid, '},\ny tick label style={\n/pgf/number format/.cd,\n');
+if o.yTickLabelFixed
+    fprintf(fid, 'fixed,\n');
+end
+if o.yTickLabelZeroFill
+    fprintf(fid, 'zerofill,\n');
+end
+fprintf(fid, 'precision=%d,\n/tikz/.cd\n},\n', o.yTickLabelPrecision);
+fprintf(fid, 'x tick label style={rotate=%f', o.xTickLabelRotation);
 if o.xTickLabelRotation ~= 0
     fprintf(fid, ',anchor=%s', o.xTickLabelAnchor);
 end
@@ -176,21 +184,8 @@ if ~isempty(o.miscTikzAxisOptions)
 end
 fprintf(fid, ']\n');
 
-if o.showZeroline
-    fprintf(fid, '%%zeroline\n\\addplot[black,line width=.5,forget plot] coordinates {(1,0)(%d,0)};\n',dd.ndat);
-end
-
-for i=1:ne
-    o.series{i}.writeSeriesForGraph(fid, dd);
-    if o.showLegend
-        le = o.series{i}.getNameForLegend();
-        if ~isempty(le)
-            fprintf(fid, '\\addlegendentry{%s}\n', le);
-        end
-    end
-end
-
 if ~isempty(o.shade)
+    fprintf(fid, '%%shading\n');
     stringsdd = strings(dd);
     x1 = find(strcmpi(date2string(o.shade(1)), stringsdd));
     x2 = find(strcmpi(date2string(o.shade(end)), stringsdd));
@@ -216,6 +211,19 @@ if ~isempty(o.shade)
     end
 end
 
+if o.showZeroline
+    fprintf(fid, '%%zeroline\n\\addplot[black,line width=.5,forget plot] coordinates {(1,0)(%d,0)};\n',dd.ndat);
+end
+
+for i=1:ne
+    o.series{i}.writeSeriesForGraph(fid, dd);
+    if o.showLegend
+        le = o.series{i}.getNameForLegend();
+        if ~isempty(le)
+            fprintf(fid, '\\addlegendentry{%s}\n', le);
+        end
+    end
+end
 fprintf(fid, '\\end{axis}\n\\end{tikzpicture}%%');
 if fclose(fid) == -1
     error('@graph.writeGraphFile: closing %s\n', o.filename);
diff --git a/matlab/reports/@report/write.m b/matlab/reports/@report/write.m
index 6d1e187d1..0970f0c9e 100644
--- a/matlab/reports/@report/write.m
+++ b/matlab/reports/@report/write.m
@@ -69,11 +69,13 @@ fprintf(fid, '\\renewcommand{\\textfraction}{0.05}\n');
 fprintf(fid, '\\renewcommand{\\topfraction}{0.8}\n');
 fprintf(fid, '\\renewcommand{\\bottomfraction}{0.8}\n');
 fprintf(fid, '\\setlength{\\parindent}{0in}\n');
+fprintf(fid, '\\setlength{\\tabcolsep}{1em}\n');
 fprintf(fid, '\\newlength\\sectionheight\n');
 fprintf(fid, '\\begin{document}\n');
 fprintf(fid, '\\pgfdeclarelayer{background0}\n');
 fprintf(fid, '\\pgfdeclarelayer{background1}\n');
 fprintf(fid, '\\pgfsetlayers{background0,background1,main}\n');
+fprintf(fid, '\\pgfplotsset{tick scale binop={\\times},\ntrim axis left}\n');
 fprintf(fid, '\\centering\n');
 
 nps = length(o.pages);
diff --git a/matlab/utilities/dates/add_periods_to_array_of_dates.m b/matlab/utilities/dates/add_periods_to_array_of_dates.m
index 7274a1a78..48a371739 100644
--- a/matlab/utilities/dates/add_periods_to_array_of_dates.m
+++ b/matlab/utilities/dates/add_periods_to_array_of_dates.m
@@ -27,12 +27,16 @@ time(:,1) = time(:,1) + fix(p/freq);
 time(:,2) = time(:,2) + rem(p,freq);
 
 id1 = find(time(:,2)>freq);
-time(id1,1) = time(id1,1) + 1;
-time(id1,2) = time(id1,2) - freq;
+if ~isempty(id1)
+    time(id1,1) = time(id1,1) + 1;
+    time(id1,2) = time(id1,2) - freq;
+end
 
 id2 = find(time(:,2)<1);
-time(id2,1) = time(id2,1) - 1;
-time(id2,2) = time(id2,2) + freq;
+if ~isempty(id2)
+    time(id2,1) = time(id2,1) - 1;
+    time(id2,2) = time(id2,2) + freq;
+end
 
 %@test:1
 %$ t(1) = dyn_assert(add_periods_to_array_of_dates([1950 1], 4, 1),[1950 2]);
diff --git a/matlab/utilities/dates/add_periods_to_date.m b/matlab/utilities/dates/add_periods_to_date.m
new file mode 100644
index 000000000..22ff62c5d
--- /dev/null
+++ b/matlab/utilities/dates/add_periods_to_date.m
@@ -0,0 +1,57 @@
+function time = add_periods_to_date(time, freq, p)  % --*-- Unitary tests --*--
+
+% Adds a p periods (p can be negative) to a date (or a set of dates) characterized by array time and frequency freq.
+    
+% Copyright (C) 2013 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+time(1) = time(1) + fix(p/freq);
+time(2) = time(2) + rem(p,freq);
+
+if time(2)>freq
+    time(1) = time(1) + 1;
+    time(2) = time(2) - freq;
+end
+
+if time(2)<1
+    time(1) = time(1) - 1;
+    time(2) = time(2) + freq;
+end
+
+%@test:1
+%$ t(1) = dyn_assert(add_periods_to_date([1950 1], 4, 1),[1950 2]);
+%$ t(2) = dyn_assert(add_periods_to_date([1950 1], 4, 2),[1950 3]);
+%$ t(3) = dyn_assert(add_periods_to_date([1950 1], 4, 3),[1950 4]);
+%$ t(4) = dyn_assert(add_periods_to_date([1950 1], 4, 4),[1951 1]);
+%$ t(5) = dyn_assert(add_periods_to_date([1950 1], 4, 5),[1951 2]);
+%$ t(6) = dyn_assert(add_periods_to_date([1950 1], 4, 6),[1951 3]);
+%$ t(7) = dyn_assert(add_periods_to_date([1950 1], 4, 7),[1951 4]);
+%$ t(8) = dyn_assert(add_periods_to_date([1950 1], 4, 8),[1952 1]);
+%$ T = all(t);
+%@eof:1
+
+%@test:2
+%$ t(1) = dyn_assert(add_periods_to_date([1950 1], 4, -1),[1949 4]);
+%$ t(2) = dyn_assert(add_periods_to_date([1950 1], 4, -2),[1949 3]);
+%$ t(3) = dyn_assert(add_periods_to_date([1950 1], 4, -3),[1949 2]);
+%$ t(4) = dyn_assert(add_periods_to_date([1950 1], 4, -4),[1949 1]);
+%$ t(5) = dyn_assert(add_periods_to_date([1950 1], 4, -5),[1948 4]);
+%$ t(6) = dyn_assert(add_periods_to_date([1950 1], 4, -6),[1948 3]);
+%$ t(7) = dyn_assert(add_periods_to_date([1950 1], 4, -7),[1948 2]);
+%$ t(8) = dyn_assert(add_periods_to_date([1950 1], 4, -8),[1948 1]);
+%$ T = all(t);
+%@eof:2
\ No newline at end of file
diff --git a/matlab/utilities/dates/isdate.m b/matlab/utilities/dates/isdate.m
index ebcd2816e..2fb603803 100644
--- a/matlab/utilities/dates/isdate.m
+++ b/matlab/utilities/dates/isdate.m
@@ -25,11 +25,11 @@ function b = isdate(str)  % --*-- Unitary tests --*--
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-if length(str)>1
+%if length(str)>1
     b = isquaterly(str) || isyearly(str) || ismonthly(str) || isweekly(str);
-else
-    b = 0;
-end
+    %else
+    %b = 0;
+    %end
 
 %@test:1
 %$
diff --git a/matlab/utilities/dates/ismonthly.m b/matlab/utilities/dates/ismonthly.m
index 7d3a60f65..a66d06e69 100644
--- a/matlab/utilities/dates/ismonthly.m
+++ b/matlab/utilities/dates/ismonthly.m
@@ -26,7 +26,11 @@ function b = ismonthly(str)  % --*-- Unitary tests --*--
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
 if ischar(str)
-    b = ~isempty(regexp(str,'^-?[0-9]*[Mm]([1-9]|1[0-2])$'));
+    if isempty(regexp(str,'^-?[0-9]*[Mm]([1-9]|1[0-2])$','once'))
+        b = 0;
+    else
+        b = 1;
+    end
 else
     b = 0;
 end
diff --git a/matlab/utilities/dates/isquaterly.m b/matlab/utilities/dates/isquaterly.m
index 92b1c6dc2..2b154ef3c 100644
--- a/matlab/utilities/dates/isquaterly.m
+++ b/matlab/utilities/dates/isquaterly.m
@@ -25,8 +25,12 @@ function b = isquaterly(str)  % --*-- Unitary tests --*--
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-if ischar(str)    
-    b = ~isempty(regexp(str,'^-?[0-9]*[Qq][1-4]$'));
+if ischar(str)
+    if isempty(regexp(str,'^-?[0-9]*[Qq][1-4]$','once'))
+        b = 0;
+    else
+        b = 1;
+    end
 else
     b = 0;
 end
diff --git a/matlab/utilities/dates/isweekly.m b/matlab/utilities/dates/isweekly.m
index ea7eca992..d9988f525 100644
--- a/matlab/utilities/dates/isweekly.m
+++ b/matlab/utilities/dates/isweekly.m
@@ -26,7 +26,11 @@ function b = isweekly(str)  % --*-- Unitary tests --*--
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
 if ischar(str)
-    b = ~isempty(regexp(str,'^-?[0-9]*[Ww]([1-9]|[1-4][0-9]|5[0-2])$'));
+    if isempty(regexp(str,'^-?[0-9]*[Ww]([1-9]|[1-4][0-9]|5[0-2])$','once'))
+        b = 0;
+    else
+        b = 1;
+    end
 else
     b = 0;
 end
diff --git a/matlab/utilities/dates/isyearly.m b/matlab/utilities/dates/isyearly.m
index ef56c46b9..d1345a0a1 100644
--- a/matlab/utilities/dates/isyearly.m
+++ b/matlab/utilities/dates/isyearly.m
@@ -26,7 +26,11 @@ function b = isyearly(str)  % --*-- Unitary tests --*--
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
 if ischar(str)
-    b = ~isempty(regexp(str,'^-?[0-9]*[YyAa]$'));
+    if isempty(regexp(str,'^-?[0-9]*[YyAa]$','once'))
+        b = 0;
+    else
+        b = 1;
+    end
 else
     b = 0;
 end
diff --git a/matlab/utilities/dates/string2date.m b/matlab/utilities/dates/string2date.m
index 139b65dbf..bb91a05fb 100644
--- a/matlab/utilities/dates/string2date.m
+++ b/matlab/utilities/dates/string2date.m
@@ -26,7 +26,7 @@ end
 if isyearly(a)
     year = 1:(regexp(a,'[AaYy]')-1);
     date.freq = 1;
-    date.time = write_time_field(a, year);
+    date.time = write_time_field_y(a, year);
     return
 end
 
@@ -51,13 +51,14 @@ if isweekly(a)
     return
 end
 
+
 function b = write_time_field(c, d)
-    b(1) = str2num(c(d));
-    if ismember(c(d(end)+1),{'Y','y','A','a'})
-        b(2) = 1;
-    else
-        b(2) = str2num(c(d(end)+2:end));
-    end
+    b(1) = str2double(c(d));
+    b(2) = str2double(c(d(end)+2:end));
+
+function b = write_time_field_y(c, d)
+    b(1) = str2double(c(d));
+    b(2) = 1;
 
 %@test:1
 %$
diff --git a/matlab/utilities/dseries/from.m b/matlab/utilities/dseries/from.m
new file mode 100644
index 000000000..9351a02e3
--- /dev/null
+++ b/matlab/utilities/dseries/from.m
@@ -0,0 +1,424 @@
+function from(varargin)
+
+% Copyright (C) 2014 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+lvarargin = lower(varargin);
+
+to_id = strmatch('to',lvarargin);
+do_id = strmatch('do',lvarargin);
+
+if isempty(to_id) || isempty(do_id)
+    error(get_error_message_0())
+end
+
+if do_id<to_id
+    msg = sprinf('dseries::from: Wrong syntax! The TO keyword must preceed the DO keyword.\n');
+    error(get_error_message_0(msg))
+end
+
+if ~isdate(varargin{1})
+    % The first argument is not a string formatted date. Test if this argument refers to a dates object
+    % in the caller workspace.
+    try
+        d1 = evalin('caller', varargin{1});
+        if ~isdates(d1)
+            error(['dseries::from: Variable ' varargin{1} ' is not a dates object!'])
+        end
+    catch
+        error(['dseries::from: Variable ' varargin{1} ' is unknown!'])
+    end
+    if ~exist('d1')
+        msg = sprintf('dseries::from: Wrong syntax! The FROM statement must be followed by a string formatted date.\n');
+        error(get_error_message_0(msg))
+    end
+else
+    d1 = dates(varargin{1}); % First date
+end
+
+if ~isequal(to_id,2)
+    msg = sprintf('dseries::from: Wrong syntax! The first dates object must be immediately followed by the TO keyword.\n');
+    error(get_error_message_0(msg))
+end
+
+if ~isdate(varargin{3})
+    % The third argument is not a string formatted date. Test if this argument refers to a dates object
+    % in the caller workspace.
+    try
+        d2 = evalin('caller', varargin{3});
+        if ~isdates(d2)
+            error(['dseries::from: Variable ' varargin{3} ' is not a dates object!'])
+        end
+    catch
+        error(['dseries::from: Variable ' varargin{3} ' is unknown!'])
+    end
+    if ~exist('d2')
+        msg = sprintf('dseries::from: Wrong syntax! The TO keyword must be followed by a second dates object.\n');
+        error(get_error_message_0(msg))
+    end
+else
+    d2 = dates(varargin{3}); % Last date
+end
+
+if d1>d2
+    error('dseries::from: The first date must preceed the second one!')
+end
+
+if ~isequal(do_id,4)
+    msg = sprintf('dseries::from: Wrong syntax! The second dates object must be immediately followed by the DO keyword.\n');
+    error(get_error_message_0(msg))
+end
+
+% Build the recursive expression.
+EXPRESSION = char([varargin{5:end}]);
+
+% Check that the expression is an assignment
+equal_id = strfind(EXPRESSION,'=');
+if isempty(equal_id)
+    error('dseries::from: Wrong syntax! The expression following the DO keyword must be an assignment (missing equal symbol).')
+end
+
+% Issue ann error message if the user attempts to do more than one assignment.
+if ~isequal(length(equal_id),1)
+    error('dseries::from: Not yet implemented! Only one assignment is allowed in the FROM-TO-DO statement.')
+end
+
+% Get all the variables involved in the recursive expression.
+variables = unique(regexpi(EXPRESSION, '\w*\(t\)|\w*\(t\-\d\)|\w*\(t\+\d\)|\w*\.\w*\(t\)|\w*\.\w*\(t\-\d\)|\w*\.\w*\(t\+\d\)','match'));
+
+% Copy EXPRESSION in expression. In the next loop we will remove all indexed variables from expression.
+expression = EXPRESSION;
+
+% Build an incidence table (max lag/lead for each variable)
+%
+% Column 1: Name of the variable.
+% Column 2: Maximum lag order.
+% Column 3: Equal to 1 if the variable appears at the current period, 0 otherwise.
+% Column 4: Maximum lead order.
+% Column 5: Vector of effective lag orders.
+% Column 6: Vector of effective lead orders.
+%
+% Initialization.
+leadlagtable = cell(0,6);
+% Loop over the variables (dseries objects).
+for i=1:length(variables)
+    expression = strrep(expression,variables{i},'');
+    current = ~isempty(regexpi(variables{i},'\(t\)'));
+    lag = ~isempty(regexpi(variables{i},'\(t\-\d\)'));
+    lead = ~isempty(regexpi(variables{i},'\(t\+\d\)'));
+    start = regexpi(variables{i},'\(t\)|\(t\-\d\)|\(t\+\d\)');
+    index = variables{i}(start:end);
+    variables(i) = {variables{i}(1:start-1)};
+    if isempty(leadlagtable)
+        leadlagtable(1,1) = {variables{i}};
+        if current
+            leadlagtable(1,3) = {1};
+        else
+            leadlagtable(1,3) = {0};
+        end
+        if lag
+            tmp = regexpi(index,'\d','match');
+            leadlagtable(1,2) = {str2double(tmp{1})};
+            leadlagtable(1,5) = {str2double(tmp{1})};
+        else
+            leadlagtable(1,2) = {0};
+            leadlagtable(1,5) = {[]};
+        end
+        if lead
+            tmp = regexpi(index,'\d','match');
+            leadlagtable(1,4) = {str2double(tmp{1})};
+            leadlagtable(1,6) = {str2double(tmp{1})};
+        else
+            leadlagtable(1,4) = {0};
+            leadlagtable(1,6) = {[]};
+        end
+    else
+        linea = strmatch(variables{i},leadlagtable(:,1));
+        if isempty(linea)
+            % This is a new variable!
+            linea = size(leadlagtable,1)+1;
+            leadlagtable(linea,1) = {variables{i}};
+            leadlagtable(linea,2) = {0};
+            leadlagtable(linea,3) = {0};
+            leadlagtable(linea,4) = {0};
+            leadlagtable(linea,5) = {[]};
+            leadlagtable(linea,6) = {[]};
+        end
+        if current
+            leadlagtable(linea,3) = {1};
+        end
+        if lag
+            tmp = regexpi(index,'\d','match');
+            leadlagtable(linea,2) = {max(str2double(tmp{1}),leadlagtable{linea,2})};
+            leadlagtable(linea,5) = {sortrows([leadlagtable{linea,5}; str2double(tmp{1})])};
+        end
+        if lead
+            tmp = regexpi(index,'\d','match');
+            leadlagtable(linea,4) = {max(str2double(tmp{1}),leadlagtable{linea,4})};
+            leadlagtable(linea,6) = {sortrows([leadlagtable{linea,6}; str2double(tmp{1})])};
+        end
+    end
+end
+
+% Set the number of variables
+number_of_variables = size(leadlagtable,1);
+
+% Initialize a cell array containing the names of the variables.
+variable_names = cell(1);
+
+% Test that all the involved variables are available dseries objects. Also check that
+% these time series are defined over the time range given by d1 and d2 (taking care of
+% the lags and leads) and check that each object is a singleton
+for i=1:number_of_variables
+    current_variable = leadlagtable{i,1};
+    idvar = strfind(current_variable,'.');
+    if isempty(idvar)
+        idvar = 0;
+    end
+    if idvar
+        current_variable_0 = current_variable(1:idvar-1);
+    else
+        current_variable_0 = current_variable;
+    end
+    try
+        var = evalin('caller',current_variable_0);
+    catch
+        error(['dseries::from: Variable ' current_variable_0 ' is unknown!'])
+    end
+    if idvar
+        try
+            eval(sprintf('var = var.%s;',current_variable(idvar+1:end)))
+        catch
+            error(sprintf('dseries::from: Variable %s is not a member of dseries oject %s!', current_variable(idvar+1:end), current_variable_0))
+        end
+    end
+    if ~isdseries(var)
+        error(['dseries::from: Variable ' current_variable ' is not a dseries object!'])
+    else
+        if ~var.vobs
+            msg = sprintf('dseries::from: Object %s must not be empty!\n',current_variable);
+            error(msg)
+        end
+        if var.vobs>1
+            msg = sprintf('dseries::from: Object %s must contain only one variable!\n',current_variable);
+            error(msg)
+        end
+        if i>1
+            if ismember(var.name,variable_names)
+                error('dseries::from: All the dseries objects should contain variables with different names!')
+            else
+                variable_names(i) = {var.name{1}};
+            end
+        else
+            variable_names(i) = {var.name{1}};
+        end
+        if d1<var.dates(1)+leadlagtable{i,2}
+            msg = sprintf('dseries::from: Initial date of the loop (%s) is inconsistent with %s''s range!\n',char(d1),current_variable);
+            msg = [msg, sprintf('               Initial date should be greater than or equal to %s.',char(var.dates(1)+leadlagtable{i,2}))];
+            error(msg)
+        end
+        if d2>var.dates(end)-leadlagtable{i,4}
+            % The first variable should be the assigned variable (will be tested later)
+            if ~isassignedvariable(leadlagtable{i,1},EXPRESSION)
+                msg = sprintf('dseries::from: Last date of the loop (%s) is inconsistent with %s''s range!\n',char(d2),current_variable);
+                msg = [msg, sprintf('               Last date should be less than or equal to %s.',char(var.dates(end)-leadlagtable{i,4}))];
+                error(msg)
+            else
+                var = [var; dseries(NaN((d2-var.dates(end)),1),var.dates(end)+1:d2,var.name)];
+            end
+        end
+        eval(sprintf('%s = var;',current_variable));
+    end
+end
+
+% Get the name of the assigned variable (with time index)
+assignedvariablename = regexpi(EXPRESSION(1:equal_id-1), '\w*\(t\)|\w*\(t\-\d\)|\w*\(t\+\d\)|\w*\.\w*\(t\)|\w*\.\w*\(t\-\d\)|\w*\.\w*\(t\+\d\)','match');
+if isempty(assignedvariablename)
+    error('dseries::from: Wrong syntax! The expression following the DO keyword must be an assignment (missing variable before the equal symbol).')
+end
+if length(assignedvariablename)>1
+    error('dseries::from: No more than one variable can be assigned!')
+end
+% Check if the model is static
+start = regexpi(assignedvariablename{1},'\(t\)|\(t\-\d\)|\(t\+\d\)');
+index = assignedvariablename{1}(start:end);
+assignedvariablename = assignedvariablename{1}(1:start-1);
+indva = strmatch(assignedvariablename, leadlagtable(:,1));
+dynamicmodel = ~isempty(regexpi(EXPRESSION(equal_id:end), ...
+                                sprintf('%s\\(t\\)|%s\\(t\\-\\d\\)|%s\\(t\\+\\d\\)',assignedvariablename,assignedvariablename,assignedvariablename),'match'));
+% Check that the dynamic model for the endogenous variable is not forward looking.
+if dynamicmodel
+    indum = index2num(index);
+    if indum<leadlagtable{indva,4}
+        error('dseries::from: It is not possible to simulate a forward looking model!')
+    end
+end
+% Check that the assigned variable does not depend on itself (the assigned variable can depend on its past level but not on the current level).
+if dynamicmodel
+    tmp = regexpi(EXPRESSION(equal_id+1:end), ...
+                  sprintf('%s\\(t\\)|%s\\(t\\-\\d\\)|%s\\(t\\+\\d\\)',assignedvariablename,assignedvariablename,assignedvariablename),'match');
+    tmp = cellfun(@extractindex, tmp);
+    tmp = cellfun(@index2num, tmp);
+    if ~all(tmp(:)<indum)
+        error(sprintf('dseries::from: On the righthand side, the endogenous variable, %s, must be indexed by %s at most.',assignedvariablename,num2index(indum-1)))
+    end
+end
+
+% Put all the variables in a unique dseries object.
+list_of_variables = leadlagtable{1,1};
+for i=2:number_of_variables
+   list_of_variables = [list_of_variables, ',' leadlagtable{i,1}];
+end
+eval(sprintf('tmp = [%s];', list_of_variables));
+
+% Get base time index
+t1 = find(d1==tmp.dates);
+t2 = find(d2==tmp.dates);
+
+% Get data
+data = tmp.data;
+
+% Isolate the (potential) parameters in the expression to be evaluated
+[~, TMP314] = strsplit(expression,'([0-9]*\.[0-9]*|\w*)','DelimiterType','RegularExpression','CollapseDelimiters',false);
+% Here I remove the numbers (TMP314 -> TMP314159).
+TMP3141 = regexp(TMP314,'(([0-9]*\.[0-9]*)|([0-9]*))','match');
+TMP31415 = find(cellfun(@isempty,TMP3141));
+TMP314159 = TMP314(TMP31415);
+
+if dynamicmodel
+    % Transform EXPRESSION by replacing calls to the dseries objects by references to data.
+    for i=1:number_of_variables
+        EXPRESSION = regexprep(EXPRESSION,sprintf('%s\\(t\\)',leadlagtable{i,1}),sprintf('data(t,%s)',num2str(i)));
+        for j=1:length(leadlagtable{i,5})
+            lag = leadlagtable{i,5}(j);
+            EXPRESSION = regexprep(EXPRESSION,sprintf('%s\\(t-%s\\)',leadlagtable{i,1},num2str(lag)), ...
+                                   sprintf('data(t-%s,%s)',num2str(lag),num2str(i)));
+        end
+        for j=1:length(leadlagtable{i,6})
+            lead = leadlagtable{i,6}(j);
+            EXPRESSION = regexprep(EXPRESSION,sprintf('%s\\(t+%s\\)',leadlagtable{i,1},num2str(lead)), ...
+                                   sprintf('data(t+%s,%s)',num2str(lead),num2str(i)));
+        end
+    end
+    % Get values for the parameters (if any)
+    if ~isempty(TMP314159)
+        for i=1:length(TMP314159)
+            wordcandidate = TMP314159{i};
+            try % If succesful, word is a reference to a variable in the caller function/script.
+                thiswordisaparameter = evalin('caller', wordcandidate);
+                eval(sprintf('%s = thiswordisaparameter;',wordcandidate));
+            catch
+                % I assume that word is a reference to a function.
+            end
+        end
+    end
+    % Do the job. Evaluate the recursion.
+    eval(sprintf('for t=%s:%s, %s; end',num2str(t1),num2str(t2),EXPRESSION));
+else
+    % Transform EXPRESSION by replacing calls to the dseries objects by references to data.
+    for i=1:number_of_variables
+        EXPRESSION = regexprep(EXPRESSION,sprintf('%s\\(t\\)',leadlagtable{i,1}), ...
+                               sprintf('data(%s:%s,%s)',num2str(t1),num2str(t2),num2str(i)));
+        for j=1:length(leadlagtable{i,5})
+            lag = leadlagtable{i,5}(j);
+            EXPRESSION = regexprep(EXPRESSION,sprintf('%s\\(t-%s\\)',leadlagtable{i,1},num2str(lag)), ...
+                                   sprintf('data(%s:%s,%s)',num2str(t1-lag),num2str(t2-lag),num2str(i)));
+        end
+        for j=1:length(leadlagtable{i,6})
+            lead = leadlagtable{i,6}(j);
+            EXPRESSION = regexprep(EXPRESSION,sprintf('%s\\(t+%s\\)',leadlagtable{i,1},num2str(lead)), ...
+                                   sprintf('data(%s:%s,%s)',num2str(t1-lead),num2str(t2-lead),num2str(i)));
+        end
+    end
+    % Transform some operators (^ -> .^, / -> ./ and * -> .*)
+    EXPRESSION = strrep(EXPRESSION,'^','.^');
+    EXPRESSION = strrep(EXPRESSION,'*','.*');
+    EXPRESSION = strrep(EXPRESSION,'/','./');
+    % Get values for the parameters (if any)
+    if ~isempty(TMP314159)
+        for i=1:length(TMP314159)
+            wordcandidate = TMP314159{i};
+            try % If succesful, word is a reference to a variable in the caller function/script.
+                thiswordisaparameter = evalin('caller', wordcandidate);
+                eval(sprintf('%s = thiswordisaparameter;',wordcandidate));
+            catch
+                % I assume that word is a reference to a function.
+            end
+        end
+    end
+    % Do the job. Evaluate the static expression.
+    eval(sprintf('%s;',EXPRESSION));
+end
+
+% Put assigned variable back in the caller workspace...
+if isempty(strfind(assignedvariablename,'.'))
+    eval(sprintf('assignin(''caller'', ''%s'', dseries(data(:,indva),%s.init,%s.name,%s.tex));', ...
+                 assignedvariablename,assignedvariablename,assignedvariablename,assignedvariablename))
+else
+    DATA = num2cell(data(:,indva));
+    strdata = sprintf('%f ', DATA{:});
+    evalin('caller',sprintf('%s = dseries(transpose([%s]),%s.init,%s.name,%s.tex);', ...
+                            assignedvariablename,strdata,assignedvariablename,assignedvariablename,assignedvariablename))
+end
+
+function msg = get_error_message_0(msg)
+    if ~nargin
+        msg = sprintf('Wrong syntax! The correct syntax is:\n\n');
+    else
+        msg = [msg, sprintf('The correct syntax is:\n\n')];
+    end
+    msg = [msg, sprintf('    from d1 to d2 do SOMETHING\n\n')];
+    msg = [msg, sprintf('where d1<d2 are dates objects, and SOMETHING is a recursive expression involving dseries objects.')];
+
+
+function index = extractindex(str)
+    index = regexpi(str,'\(t\)|\(t\-\d\)|\(t\+\d\)','match');
+
+
+function i = index2num(id)
+    if isequal('(t)',id)
+        i = 0;
+        return
+    end
+    if isequal('-',id(3))
+        i = - str2num(id(4:end-1));
+    else
+        i = str2num(id(4:end-1));
+    end
+
+
+function id = num2index(i)
+    if isequal(i,0)
+        id = '(t)';
+        return
+    end
+    if i<0
+        id = ['(t-' int2str(abs(i)) ')'];
+    else
+        id = ['(t+' int2str(i) ')'];
+    end
+
+function i = isassignedvariable(var,expr)
+    idv = strfind(expr,var);
+    idq = strfind(expr,'=');
+    if ~isempty(idv)
+        if idv(1)<idq
+            i = 1;
+            return
+        end
+    end
+    i = 0;
\ No newline at end of file
diff --git a/matlab/utilities/general/name2tex.m b/matlab/utilities/general/name2tex.m
index 61ec30de0..0bf02a2cb 100644
--- a/matlab/utilities/general/name2tex.m
+++ b/matlab/utilities/general/name2tex.m
@@ -1,7 +1,7 @@
 function tex = name2tex(name, info) % --*-- Unitary tests --*--
 
 % Converts plain text name into tex name.
- 
+
 % Copyright (C) 2013 Dynare Team
 %
 % This file is part of Dynare.
@@ -64,7 +64,7 @@ if info
             if gotonextcondition && ndx>2
                 if idx(end)<ntx
                     texname = [ texname(1:idx(end)-2) '_{' texname(idx(end)+1:end) '}' ];
-                else 
+                else
                     texname = [ texname(1:idx(end-1)-2) '_{' texname(idx(end-1)+1:end) '}' ];
                 end
             end
@@ -93,9 +93,9 @@ end
 %$ t14 = name2tex('azert_uiop',1);
 %$ t15 = name2tex('azert_uiop_qsdfg',1);
 %$ t16 = name2tex('azert_uiop_qsdfg_',1);
-%$ t17 = name2tex('_azert_uiop_qsdfg',1); 
+%$ t17 = name2tex('_azert_uiop_qsdfg',1);
 %$ t18 = name2tex('_azert_uiop_qsdfg_',1);
-%$ 
+%$
 %$ t(1) = dyn_assert(strcmp(t1,'\\_azert'),1);
 %$ t(2) = dyn_assert(strcmp(t2,'azert\\_'),1);
 %$ t(3) = dyn_assert(strcmp(t3,'\\_azert\\_'),1);
@@ -106,7 +106,7 @@ end
 %$ t(8) = dyn_assert(strcmp(t8,'\\_azert\\_uiop\\_qsdfg\\_'),1);
 %$ t(9) = dyn_assert(strcmp(t11,'\\_azert'),1);
 %$ t(10) = dyn_assert(strcmp(t12,'azert\\_'),1);
-%$ t(11) = dyn_assert(strcmp(t13,'\\_azert\\_'),1); 
+%$ t(11) = dyn_assert(strcmp(t13,'\\_azert\\_'),1);
 %$ t(12) = dyn_assert(strcmp(t14,'azert_{uiop}'),1);
 %$ t(13) = dyn_assert(strcmp(t15,'azert\\_uiop_{qsdfg}'),1);
 %$ t(14) = dyn_assert(strcmp(t16,'azert\\_uiop_{qsdfg\\_}'),1);
@@ -185,7 +185,6 @@ end
 %@eof:3
 
 %@test:4
-%$ pwd
 %$ try
 %$     db = dseries('csv/dd.csv');
 %$     t(1) = 1;
diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc
index 369f46b93..891612146 100644
--- a/preprocessor/ModFile.cc
+++ b/preprocessor/ModFile.cc
@@ -554,6 +554,10 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool no_log, b
     mOutputFile << "clear all" << endl;
 
   mOutputFile << "tic;" << endl
+	      << "% Save empty dates and dseries objects in memory." << endl
+	      << "dates('initialize');" << endl
+	      << "dseries('initialize');" << endl
+	      << "% Define global variables." << endl
               << "global M_ oo_ options_ ys0_ ex0_ estimation_info" << endl
               << "options_ = [];" << endl
               << "M_.fname = '" << basename << "';" << endl
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 2600f1295..2e8b02fa9 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -177,7 +177,8 @@ MODFILES = \
 	loglinear/example4_loglinear.mod \
 	smoother2histval/fs2000_simul.mod \
 	smoother2histval/fs2000_smooth.mod \
-	smoother2histval/fs2000_smooth_stoch_simul.mod
+	smoother2histval/fs2000_smooth_stoch_simul.mod \
+	reporting/example1.mod
 
 XFAIL_MODFILES = ramst_xfail.mod \
 	estim_param_in_shock_value.mod
diff --git a/tests/reporting/example1.mod b/tests/reporting/example1.mod
new file mode 100644
index 000000000..bec34f5b5
--- /dev/null
+++ b/tests/reporting/example1.mod
@@ -0,0 +1,81 @@
+// Example 1 from Collard's guide to Dynare
+var y, c, k, a, h, b;
+varexo e, u;
+
+verbatim;
+% I want these comments included in
+% example1.m
+%
+var = 1;
+end;
+
+parameters beta, rho, alpha, delta, theta, psi, tau;
+
+alpha = 0.36;
+rho   = 0.95;
+tau   = 0.025;
+beta  = 0.99;
+delta = 0.025;
+psi   = 0;
+theta = 2.95;
+
+phi   = 0.1;
+
+model;
+c*theta*h^(1+psi)=(1-alpha)*y;
+k = beta*(((exp(b)*c)/(exp(b(+1))*c(+1)))
+    *(exp(b(+1))*alpha*y(+1)+(1-delta)*k));
+y = exp(a)*(k(-1)^alpha)*(h^(1-alpha));
+k = exp(b)*(y-c)+(1-delta)*k(-1);
+a = rho*a(-1)+tau*b(-1) + e;
+b = tau*a(-1)+rho*b(-1) + u;
+end;
+
+initval;
+y = 1.08068253095672;
+c = 0.80359242014163;
+h = 0.29175631001732;
+k = 11.08360443260358;
+a = 0;
+b = 0;
+e = 0;
+u = 0;
+end;
+
+shocks;
+var e; stderr 0.009;
+var u; stderr 0.009;
+var e, u = phi*0.009*0.009;
+end;
+
+stoch_simul;
+
+shocke = dseries();
+shocku = dseries();
+@#define endovars=["y", "c", "k", "a", "h", "b"]
+@#for var in endovars
+  shocke = [shocke dseries(@{var}_e, 2014q3, '@{var}')];
+  shocku = [shocku dseries(@{var}_u, 2014q3, '@{var}')];
+@#endfor
+
+r = report();
+@#for shock in ["e", "u"]
+    r = r.addPage('title',{'Dseries/Report Example','Shock @{shock}'},...
+                  'titleFormat', {'\Large\bfseries', '\large\bfseries'});
+    r = r.addSection('cols', 2);
+@#  for var in endovars
+      r = r.addGraph('data', shock@{shock}.@{var}, 'title', '@{var}', ...
+                     'showGrid', false, 'yTickLabelPrecision', 2, 'yTickLabelZeroFill', false);
+      r = r.addSeries('graphHline', 0, 'graphLineColor', 'red');
+@#  endfor
+    r = r.addVspace('number', 2);
+    r = r.addSection('cols', 1);
+    r = r.addTable('range', 2022q1:2024q1, 'precision', 5);
+
+@#  for var in endovars
+      r = r.addSeries('data', shock@{shock}.@{var});
+@#  endfor
+@#endfor
+
+r.write();
+r.compile();
\ No newline at end of file