LTSts

LTSts extends LTS estimator to time series

Syntax

Description

It is possible to set a model with a trend (up to third order), a seasonality (constant or of varying amplitude and with a different number of harmonics) and a level shift (in this last case it is possible to specify the window in which level shift has to be searched for).

example

out =LTSts(y) Simulated data with linear trend and level shift.

example

out =LTSts(y, Name, Value) Airline data: linear trend + just one harmonic for seasonal component.

example

[out, varargout] =LTSts(___) Model with linear trend and six harmonics for seasonal component.

Examples

expand all

  • Simulated data with linear trend and level shift.
  • No seasonal component.

        rng('default')
        n=45;
        a=1;
        b=0.8;
        sig=1;
        seq=(1:n)';
        y=a+b*seq+sig*randn(n,1);
        % Add a level shift in the simulated series
        y(round(n/2):end)=y(round(n/2):end)+10;
        % model with a linear trend, non seasonal and level shift
        model=struct;
        model.trend=1;
        model.seasonal=0;
        % Potential level shift position is investigated in positions:
        % t=10, t=11, ..., t=T-10.
        model.lshift=10;
        out=LTSts(y,'model',model,'plots',1);
        % Using the notation of the paper RPRH: A=1, B=1, G=0 and $\delta_1>0$.
        str=strcat('A=1, B=0, G=0, $\delta_2=',num2str(out.posLS),'$');
        title(findobj(gcf,'-regexp','Tag','LTSts:ts'),str,'Interpreter','Latex');
    

  • Airline data: linear trend + just one harmonic for seasonal component.
  • Load airline data.

        %   1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960
        y = [112  115  145  171  196  204  242  284  315  340  360  417    % Jan
             118  126  150  180  196  188  233  277  301  318  342  391    % Feb
             132  141  178  193  236  235  267  317  356  362  406  419    % Mar
             129  135  163  181  235  227  269  313  348  348  396  461    % Apr
             121  125  172  183  229  234  270  318  355  363  420  472    % May
             135  149  178  218  243  264  315  374  422  435  472  535    % Jun
             148  170  199  230  264  302  364  413  465  491  548  622    % Jul
             148  170  199  242  272  293  347  405  467  505  559  606    % Aug
             136  158  184  209  237  259  312  355  404  404  463  508    % Sep
             119  133  162  191  211  229  274  306  347  359  407  461    % Oct
             104  114  146  172  180  203  237  271  305  310  362  390    % Nov
             118  140  166  194  201  229  278  306  336  337  405  432 ]; % Dec
        % Source:
        % http://datamarket.com/data/list/?q=provider:tsdl
    
        y=(y(:));
        yr = repmat((1949:1960),12,1);
        mo = repmat((1:12)',1,12);
        time = datestr(datenum(yr(:),mo(:),1));
        ts = timeseries(y(:),time,'name','AirlinePassengers');
        ts.TimeInfo.Format = 'dd-mmm-yyyy';
        tscol = tscollection(ts);
        % plot airline data
        plot(ts)
        % linear trend + just one harmonic for seasonal component
        model=struct;
        model.trend=1;              % linear trend
        model.s=12;                 % monthly time series
        model.seasonal=1;           % just one harmonic
        model.lshift=0;             % no level shift
        out=LTSts(y,'model',model,'dispresults',true);
    
        close all
        % Plot real and fitted values
        plot(y);
        hold('on')
        plot(out.yhat,'red')
        legend('real values','fitted values','Location','SouthEast')
        numpar = {'model parameters:' , 'A=1, B=1, G=0, $\delta_1=0$'};
        title(gca,numpar,'Interpreter','Latex');
    

  • Model with linear trend and six harmonics for seasonal component.
  •     model=struct;
        model.trend=1;              % linear trend
        model.s=12;                 % monthly time series
        model.seasonal=6;           % six harmonics
        model.lshift=0;             % no level shift
        out=LTSts(y,'model',model);
    
        close all
        % Plot real and fitted values
        plot(y);
        hold('on')
        plot(out.yhat,'red')
        legend('real values','fitted values','Location','SouthEast')
        numpar = {'model parameters:' , 'A=1, B=6, G=0, $\delta_1=0$'};
        title(gca,numpar,'Interpreter','Latex');
    
    

    Related Examples

  • Model with linear trend, two harmonics for seasonal component and varying amplitude using a linear trend.
  •     model=struct;
        model.trend=1;              % linear trend
        model.s=12;                 % monthly time series
        model.seasonal=102;         % two harmonics with time varying seasonality
        model.lshift=0;             % no level shift
        out=LTSts(y,'model',model);
    
        close all
        % Plot real and fitted values
        plot(y);
        hold('on')
        plot(out.yhat,'red')
        legend('real values','fitted values','Location','SouthEast')
        numpar = {'model parameters:' , 'A=1, B=2, G=1, $\delta_1=0$'};
       title(gca,numpar,'Interpreter','Latex');
    

  • Model with linear trend, six harmonics for seasonal component and varying amplitude using a linear trend).
  •     model=struct;
        model.trend=1;              % linear trend
        model.s=12;                 % monthly time series
        model.seasonal=106;         % six harmonics with linear time varying seasonality
        model.lshift=0;             % no level shift
        % out=fitTSLS(y,'model',model);
        out=LTSts(y,'model',model);
        close all
        % Plot real and fitted values
        plot(y);
        hold('on')
        plot(out.yhat,'red')
        legend('real values','fitted values','Location','SouthEast')
        numpar = {'model parameters:' , 'A=1, B=6, G=1, $\delta_1=0$'};
       title(gca,numpar,'Interpreter','Latex');
    
    

  • Contaminated time series with upward level shift.
  • Model with linear trend, six harmonics for seasonal component and varying amplitude using a linear trend).

        yLS=y;
        yLS(55:end)=yLS(55:end)+130;
        model=struct;
        model.trend=1;              % linear trend
        model.s=12;                 % monthly time series
        model.seasonal=1;
        model.lshift=13;            % impose level shift
        out=LTSts(yLS,'model',model);
        close all
        % Plot real and fitted values
        plot(yLS);
        hold('on')
        plot(out.yhat,'red')
        legend('real values','fitted values','Location','SouthEast')
        % Using the notation of the paper RPRH: A=1, B=1, G=0 and $\delta_1>0$.
        str=strcat('A=1, B=1, G=0, $\delta_2=',num2str(out.posLS),'$');
        numpar = {'model parameters:' , str};
        title(gca,numpar,'Interpreter','Latex');
    
    

  • Contaminated time series with downward level shift.
  • Model with linear trend, six harmonics for seasonal component and varying amplitude using a linear trend).

        yLS=y;
        yLS(35:end)=yLS(35:end)-300;
        model=struct;
        model.trend=1;              % linear trend
        model.s=12;                 % monthly time series
        model.seasonal=106;
        model.lshift=13;
        out=LTSts(yLS,'model',model);
        close all
        % Plot real and fitted values
        plot(yLS);
        hold('on')
        plot(out.yhat,'red')
        legend('real values','fitted values','Location','SouthEast')
        % Using the notation of the paper RPRH: A=1, B=6, G=1 and $\delta_1>0$.
        str=strcat('A=1, B=6, G=1, $\delta_2=',num2str(out.posLS),'$');
        numpar = {'model parameters:' , str};
        title(gca,numpar,'Interpreter','Latex');
    

  • Model with an explanatory variable using log-transformed series.
  •     y1=log(y);
        % Model with linear trend, six harmonics for seasonal component and
        % varying amplitude using a linear trend).
        model=struct;
        model.trend=1;              % linear trend
        model.s=12;                 % monthly time series
        model.seasonal=106;
        model.lshift=0;
        model.X=randn(length(y),1);
        out=LTSts(y1,'model',model);
        close all
        % Plot real and fitted values
        plot(y1);
        hold('on')
        plot(out.yhat,'red')
        legend('real values','fitted values','Location','SouthEast')
        % Using the notation of the paper RPRH: A=1, B=6, G=1 and $\delta_1>0$.
        str=strcat('A=1, B=6, G=1, $\delta_1=0$');
        numpar = {'model parameters:' , str};
        title(gca,numpar,'Interpreter','Latex');
    

  • Example 1 used in the paper RPRH.
  • Load airline data.

        %   1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960
        y = [112  115  145  171  196  204  242  284  315  340  360  417    % Jan
             118  126  150  180  196  188  233  277  301  318  342  391    % Feb
             132  141  178  193  236  235  267  317  356  362  406  419    % Mar
             129  135  163  181  235  227  269  313  348  348  396  461    % Apr
             121  125  172  183  229  234  270  318  355  363  420  472    % May
             135  149  178  218  243  264  315  374  422  435  472  535    % Jun
             148  170  199  230  264  302  364  413  465  491  548  622    % Jul
             148  170  199  242  272  293  347  405  467  505  559  606    % Aug
             136  158  184  209  237  259  312  355  404  404  463  508    % Sep
             119  133  162  191  211  229  274  306  347  359  407  461    % Oct
             104  114  146  172  180  203  237  271  305  310  362  390    % Nov
             118  140  166  194  201  229  278  306  336  337  405  432 ]; % Dec
        % Two short level shifts in opposite directions and an isolated outlier.
        % Add a level shift contamination plus some outliers.
        y1=y(:);
        y1(50:55)=y1(50:55)-300;
        y1(70:75)=y1(70:75)+300;
        y1(90:90)=y1(90:90)+300;
        % Create structure specifying model
        model=struct;
        model.trend=2;              % quadratic trend
        model.s=12;                 % monthly time series
        model.seasonal=204;         % number of harmonics
        model.lshift=40;            % position where to start monitoring level shift
        model.X='';
        % Create structure lts specifying lts options
        lshiftlocref=struct;
        % Set window length for local refinement.
        lshiftlocref.wlength=10;
        % Set tuning constant to use insde Huber rho function
        lshiftlocref.huberc=1.5;
        % Estimate the parameters
        [out]=LTSts(y1,'model',model,'nsamp',500,...
           'plots',1,'lshiftlocref',lshiftlocref,'msg',0);
        % Using the notation of the paper RPRH: A=2, B=4, G=2 and $\delta_1>0$.
        str=strcat('A=2, B=4, G=2, $\delta_2=',num2str(out.posLS),'$');
        numpar = {'model parameters:' , str};
        axeslast=findobj('-regexp','Tag','LTSts:ts');
        title(axeslast(end),numpar,'Interpreter','Latex');
    
        % generate the wedgeplot
        % wedgeplot(out,'transpose',true,'extradata',[y1 out.yhat]);
    

  • Example 2 used in the paper RPRH.
  • A persisting level shift and three isolated outliers, two of which in proximity of the level shift.

        close all
        y1=y(:);
        y1(68:end)=y1(68:end)+1300;
        y1(67)=y1(67)-600;
        y1(45)=y1(45)-800;
        y1(68:69)=y1(68:69)+800;
        % Create structure specifying model
        model=struct;
        model.trend=2;              % quadratic trend
        model.s=12;                 % monthly time series
        model.seasonal=204;         % number of harmonics
        model.lshift=40;            % position where to start monitoring level shift
        model.X='';
        % Create structure lts specifying lts options
        lshiftlocref=struct;
        % Set window length for local refinement.
        lshiftlocref.wlength=10;
        % Set tuning constant to use insde Huber rho function
        lshiftlocref.huberc=1.5;
        % Estimate the parameters
        [out, varargout]=LTSts(y1,'model',model,'nsamp',500,...
           'plots',1,'lshiftlocref',lshiftlocref,'msg',0);
    
        % Using the notation of the paper RPRH: A=2, B=4, G=2 and $\delta_1>0$.
        str=strcat('A=2, B=4, G=2, $\delta_2=',num2str(out.posLS),'$');
        numpar = {'model parameters:' , str};
        title(findobj('-regexp','Tag','LTSts:ts'),numpar,'Interpreter','Latex');
    
        % generate the wedgeplot
        % wedgeplot(out,'transpose',true,'extradata',[y1 out.yhat]);
    
    

  • Example 3 used in the paper RPRH.
  • A persisting level shift preceded and followed in the proximity by other two short level shifts, and an isolated outlier.

        y1=y(:);
        y1(50:55)=y1(50:55)-300;
        y1(68:end)=y1(68:end)-700;
        y1(70:75)=y1(70:75)+300;
        y1(90:90)=y1(90:90)+300;
        % Create structure specifying model
        model=struct;
        model.trend=2;              % quadratic trend
        model.s=12;                 % monthly time series
        model.seasonal=204;         % number of harmonics
        model.lshift=40;            % position where to start monitoring level shift
        model.X='';
        % Create structure lts specifying lts options
        lshiftlocref=struct;
        % Set window length for local refinement.
        lshiftlocref.wlength=10;
        % Set tuning constant to use insde Huber rho function
        lshiftlocref.huberc=1.5;
        close all
        % Estimate the parameters
        [out, varargout]=LTSts(y1,'model',model,'nsamp',500,...
           'plots',2,'lshiftlocref',lshiftlocref,'msg',0);
        % Using the notation of the paper RPRH: A=2, B=4, G=2 and $\delta_1>0$.
        str=strcat('A=2, B=4, G=2, $\delta_2=',num2str(out.posLS),'$');
        numpar = {'model parameters:' , str};
        title(findobj('-regexp','Tag','LTSts:ts'),numpar,'Interpreter','Latex');
    
        % generate the wedgeplot
        % wedgeplot(out,'transpose',true,'extradata',[y1 out.yhat]);
    
    
    IdleTimeout has been reached.
    Parallel pool using the 'local' profile is shutting down.
    

    Input Arguments

    expand all

    y — Time series to analyze. Vector.

    A row or a column vector with T elements, which contains the time series.

    Data Types: single| double

    Name-Value Pair Arguments

    Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside single quotes (' '). You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

    Example: 'model', model , 'intercept',1 , 'h',round(n*0,75) , 'bdp',0.4 , 'lts',lts , 'nsamp',500 , 'reftolALS',1e-05 , 'refstepsALS',20 , 'conflev',0.99 , 'plots',1 , 'SmallSampleCor',3 , 'msg',1 , 'nocheck',1 , 'lshiftlocref',lshiftlocref.typeres=2 , 'nbestindexes',5 , 'dispresults',true , 'yxsave',1

    model —model type.structure.

    A structure which specifies the model which will be used. The model structure contains the following fields:

    Value Description
    s

    scalar (length of seasonal period). For monthly data s=12 (default), for quartely data s=4, ...

    trend

    scalar (order of the trend component).

    trend = 1 implies linear trend with intercept (default), trend = 2 implies quadratic trend ...

    Admissible values for trend are, 0, 1, 2 and 3.

    In the paper RPRH to denote the order of the trend symbol A is used.

    seasonal

    scalar (integer specifying number of frequencies, i.e. harmonics, in the seasonal component. Possible values for seasonal are $1, 2, ..., [s/2]$, where $[s/2]=floor(s/2)$.

    For example:

    if seasonal =1 (default) we have:

    $\beta_1 \cos( 2 \pi t/s) + \beta_2 sin ( 2 \pi t/s)$;

    if seasonal =2 we have:

    $\beta_1 \cos( 2 \pi t/s) + \beta_2 \sin ( 2 \pi t/s) + \beta_3 \cos(4 \pi t/s) + \beta_4 \sin (4 \pi t/s)$.

    Note that when $s$ is even the sine term disappears for $j=s/2$ and so the maximum number of trigonometric parameters is $s-1$.

    If seasonal is a number greater than 100 then it is possible to specify how the seasonal component grows over time.

    For example, seasonal =101 implies a seasonal component which just uses one frequency which grows linearly over time as follows:

    $(1+\beta_3 t)\times ( \beta_1 cos( 2 \pi t/s) + \beta_2 \sin ( 2 \pi t/s))$.

    For example, seasonal =201 implies a seasonal component which just uses one frequency which grows in a quadratic way over time as follows:

    $(1+\beta_3 t + \beta_4 t^2)\times( \beta_1 \cos( 2 \pi t/s) + \beta_2 \sin ( 2 \pi t/s))$.

    seasonal =0 implies a non seasonal model.

    In the paper RPRH to denote the number of frequencies of the seasonal component symbol B is used, while symbol G is used to denote the order of the trend of the seasonal component.

    Therefore, for example, model.seasonal=201 corresponds to B=1 and G=2, while model.seasonal=3 corresponds to B=3 and G=0;

    X

    matrix of size T-by-nexpl containing the values of nexpl extra covariates which are likely to affect y.

    lshift

    scalar greater or equal than 0 which specifies whether it is necessary to include a level shift component. lshift = 0 (default) implies no level shift component. If lshift is an interger greater then 0 then it is possible to specify the moment to start considering level shifts. For example if lshift =13 then the following additional parameters are estimated $\beta_{LS1}* I(t \geq beta_{LS2})$ where $\beta_{LS1}$ is a real number and $\beta_{LS2}$ is an integer which assumes values 14, 14, ..., T-13.

    In general, the level shift which are considered are referred to times (lshift+1):(T-lshift).

    In the paper RPRH $\beta_{LS1}$ is denoted with symbol $\delta_1$, while, $\beta_{LS2}$ is denoted with symbol $\delta_2$.

    Remark: the default model is for monthly data with a linear trend (2 parameters) + seasonal component with just one harmonic (2 parameters), no additional explanatory variables and no level shift that is model=struct;

    model.s=12;

    model.trend=1;

    model.seasonal=1;

    model.X='';

    model.lshift=0;

    Using the notation of the paper RPRH we have A=1, B=1; and $\delta_1=0$.

    Example: 'model', model

    Data Types: struct

    intercept —Indicator for constant term.scalar.

    If 1, a model with constant term will be fitted (default), else no constant term will be included.

    Example: 'intercept',1

    Data Types: double

    h —The number of observations that determined the least trimmed squares estimator.scalar.

    h is an integer greater than p (number of columns of matrix X including the intercept but smaller then n. If the purpose is outlier detection than h does not have to be smaller than [0.5*(T+p+1)]. The default value of h is [0.75*T]. Note that if h is supplied input argument bdp is ignored.

    Example: 'h',round(n*0,75)

    Data Types: double

    bdp —breakdown point.scalar.

    It measures the fraction of outliers the algorithm should resist. In this case any value greater than 0 but smaller or equal than 0.5 will do fine. Please specify h or bdp, but not both.

    Example: 'bdp',0.4

    Data Types: double

    lts —structure which controls a set of options of the maximization procedure.structure.

    Structure with the following fields:

    Value Description
    refsteps

    scalar defining number of concentration steps (default = 2). refsteps = 0 means "raw-subsampling" without iterations.

    reftol

    scalar. Default value of tolerance for the refining steps The default value is 1e-6;

    bestr

    scalar defining number of "best betas" to remember from the subsamples. These will be later iterated until convergence.

    The default is 20 (10 of them are the best from previous iteration in case a level shift is present).

    refstepsbestr

    scalar defining maximum number of refining steps for each best subset (default=50).

    reftolbestr

    scalar. Default value of tolerance for the refining steps for each of the best subsets The default value is 1e-8.

    Remark: if lts is an empty value all default values of structure lts will be used.

    Example: 'lts',lts

    Data Types: struct

    nsamp —number of subsamples to extract.scalar | vector of length 2.

    Vector of length 1 or 2 which controls the number of subsamples which will be extracted to find the robust estimator. If lshift>0 then nsamp(1) controls the number of subsets which have to be extracted to find the solution for t=lshift. nsamp(2) controls the number of subsets which have to be extracted to find the solution for t=lshift+1, lshift+2, ..., T-lshift.

    Note that nsamp(2) is generally smaller than nsamp(1) because in order to compute the best solution for t=lshift+1, lshift+2, ..., T-lshift, we use the lts.bestr/2 best solutions from previous t (after shifting by one the position of the level shift in the estimator of beta). If lshift is >0 the default value of nsamp is (500 250). If lshift is >0 and nsamp is supplied as a scalar the default is to extract [nsamp/2] subsamples for t=lshift+1, lshift+2, ... Therefore, for example, in order to extract 600 subsamples for t=lshift and 300 subsamples for t= lshift+1 ... you can use nsamp =600 or nsamp=[600 300].

    The default value of nsamp is 1000;

    Remark: if nsamp=0 all subsets will be extracted.

    They will be (n choose p).

    Example: 'nsamp',500

    Data Types: double

    reftolALS —Tolerance inside ALS.scalar.

    Tolerance value of tolerance for the refining steps inside ALS routine. The default value is 1e-03.

    Example: 'reftolALS',1e-05

    Data Types: double

    refstepsALS —Maximum iterations inside ALS.scalar.

    Maximum number of iterations inside ALS routine. Default value of tolerance for the refining steps inside ALS routine. The default value is 50.

    Example: 'refstepsALS',20

    Data Types: double

    conflev —Confidence level.scalar.

    Scalar between 0 and 1 containing Confidence level which is used to declare units as outliers. Usually conflev=0.95, 0.975 0.99 (individual alpha) or 1-0.05/n, 1-0.025/n, 1-0.01/n (simultaneous alpha). Default value is 0.975.

    Example: 'conflev',0.99

    Data Types: double

    plots —Plots on the screen.scalar.

    If plots = 1, a two panel plot will be shown on the screen.

    The upper panel contains the orginal time series with fitted values. The bottom panel will contain the plot of robust residuals against index number. The confidence level which is used to draw the horizontal lines associated with the bands for the residuals is specified in input option conflev. If conflev is missing a nominal 0.975 confidence interval will be used. If plots =2 the following additional plots will be shown on the screen.

    1) Boxplot of the distribution of the lts.bestr values of the target function for each tentative level shift position;

    2) A two panel plot which shows the values of the local sum of squares varying the position of the level shift around the first tentative position keeping all the other parameters fixed. Top panel refers to Huberized residuals sum of squares and bottom panel refers to residual sum of squares.

    3) A plot which shows the indexes of the best nbestindexes solutions for each tentative level shift position.

    4) A plot which shows the relative frequency of inclusion of each unit in the best h-subset after lts.refsteps refining steps.

    5) A plot which shows the relative frequency of inclusion of each unit inside the best nbestindexes subsets which are brought to full convergence.

    The default value of plot is 0 i.e. no plot is shown on the screen.

    Example: 'plots',1

    Data Types: double

    SmallSampleCor —Small sample correction factor to control empirical size of the test.scalar equal to 1 | 2 (default) | 3 | 4.

    - If SmallSampleCor=1 in the reweighting step the nominal threshold based on $\chi^2_{0.99}$ is multiplied by the small sample correction factor which guarrantees that the empirical size of the test is equal to the nominal size.

    Given that the correction factors were obtained through simulation for a linear model, the number of explanatory which is used to compute the correction factor refers to all explanatory variables except the non linear components in the seasonal part of the model. For example, in a model with linear trend 4 seasonal harmonics + level shift and second order trend in the seasonal component the number of explanatory variables used is 11 = total number of variables -2 = 2 (linear trend) + 8 (4 seasonal harmonics) +1 (level shift).

    - If SmallSampleCor =2 Gervini and Yohai procedure is called with 'iterating' false and 'alpha' 0.99 is invoked, that is:

    weights=GYfilt(stdres,'iterating',false,'alpha',0.99);

    - If SmallSampleCor =3 Gervini and Yohai procedure called with 'iterating' true and 'alpha' 0.99 is invoked, that is:

    weights=GYfilt(stdres,'iterating',true,'alpha',0.99);

    - If SmallSampleCor =4 $\chi^2_{0.99}$ threshold is used that is:

    weights = abs(stdres)<=sqrt(chi2inv(0.99,1));

    Example: 'SmallSampleCor',3

    Data Types: double

    msg —Messages on the screen.scalar.

    Scalar which controls whether to display or not messages on the screen If msg==1 (default) messages are displayed on the screen about estimated time to compute the estimator and the warnings about 'MATLAB:rankDeficientMatrix', 'MATLAB:singularMatrix' and 'MATLAB:nearlySingularMatrix' are set to off else no message is displayed on the screen

    Example: 'msg',1

    Data Types: double

    nocheck —Check input arguments.scalar.

    If nocheck is equal to 1 no check is performed on matrix y and matrix X. Notice that y and X are left unchanged. In other words the additioanl column of ones for the intercept is not added. As default nocheck=0. The controls on h, bdp and nsamp still remain.

    Example: 'nocheck',1

    Data Types: double

    lshiftlocref —Parameters for local shift refinement.structure.

    This option is used just if model.lshift is greater then 0.

    In order to precisely identify level shift position it is necessary to consider a local sum of squares varying the position of the level shift around the first tentative position keeping all the other parameters fixed. This structure contains the following fields:

    Value Description
    wlength

    scalar greater than 0 which identifies the length of the window. The default value is 15, that is the tentative level shift position varies from tl-15, tl-15, ..., tl+14, tl+15, where tl is the best preliminary tentative level shift position.

    typeres

    scalar which identifies the type of residuals to consider. If typerres =1, the local residuals sum of squares is based on huberized (scaled) residuals (this is the default choice) else raw residuals are used.

    huberc

    tuning constant for Huber estimator just in case lshiftlocref.typeres=1. The default value is 2.

    Example: 'lshiftlocref',lshiftlocref.typeres=2

    Data Types: struct

    nbestindexes —position of the best solutions.positive integer.

    For each tentative level shift solution, it is interesenting to understand whether best solutions of target function come from subsets associated with current level shift solution or from best solutions from previous tentative level shift position. The indexes from 1 to lts.bestr/2 are associated with subsets just extracted. The indexes from lts.bestr/2+1 to lts.bestr are associated with best solutions from previous tentative level shift. More precisely:

    index lts.bestr/2+1 is associated with best solution from previous tentative level shift;

    index lts.bestr/2+2 is associated with second best solution from previous tentative level shift;

    ...

    nbestindexes is an integer which specifies how many indexes we want to store. The default value of nbestindexes is 3.

    Example: 'nbestindexes',5

    Data Types: double

    dispresults —Display results of final fit.boolean.

    If dispresults is true, labels of coefficients, estimated coefficients, standard errors, tstat and p-values are shown on the screen in a fully formatted way. The default value of dispresults is false.

    Example: 'dispresults',true

    Data Types: logical

    yxsave —store X and y.scalar.

    Scalar that is set to 1 to request that the response vector y and data matrix X are saved into the output structure out. Default is 0, i.e. no saving is done.

    Example: 'yxsave',1

    Data Types: double

    Output Arguments

    expand all

    out — description Structure

    A structure containing the following fields

    Value Description
    B

    Matrix containing estimated beta coefficients, (including the intercept when options.intercept=1) standard errors, t-stat and p-values The content of matrix B is as follows:

    1st col = beta coefficients The order of the beta coefficients is as follows:

    1) trend elements (if present). If the trend is of order two there are r+1 coefficients if the intercept is present otherwise there are just r components;

    2) linear part of seasonal component 2, 4, 6, ..., s-2, s-1 coefficients (if present);

    3) coefficients associated with the matrix of explanatory variables which have a potential effect on the time series under study (X);

    4) non linear part of seasonal component, that is varying amplitude. If varying amplitude is of order k there are k coefficients (if present);

    5) level shift component (if present). In this case there are two coefficients, the second (which is also the last element of vector beta) is an integer which specifies the time in which level shift takes place and the first (which is also the penultime element of vector beta) is a real number which identifies the magnitude of the upward (downward) level shift;

    2nd col = standard errors;

    3rd col = t-statistics;

    4th col = p values.

    h

    The number of observations that have determined the initial LTS estimator, i.e. the value of h.

    bs

    Vector containing the units with the smallest p+k squared residuals before the reweighting step, where p is the total number of the parameters in the model and p+k is smallest number of units such that the design matrix is full rank.

    out.bs can be used to initialize the forward search.

    Hsubset

    matrix of size T-by-(T-2*lshift) containing units forming best H subset for each tentative level shift which is considered.

    Units belonging to subset are given with their row number, units not belonging to subset have missing values ( Remark: T-2*lshift = length((lshift+1):(T-lshift)) ) This output is present just if input option model.lshift>0.

    posLS

    scalar associated with best tentative level shift position.

    This output is present just if input option model.lshift>0.

    numscale2

    matrix of size lts.bestr-by-(T-2*lshift) containing (in the columns the values of the lts.bestr smallest values of the target function. Target function = truncated residuals sum of squares.

    BestIndexes

    matrix of size nbestindexes-by-(T-2*lshift) containing in each column the indexes associated with the best nbestindexes solutions.

    The indexes from lts.bestr/2+1 to lts.bestr are associated with best solutions from previous tentative level shift.

    More precisely:

    index lts.bestr/2+1 is associated with best solution from previous tentative level shift;

    index lts.bestr/2+2 is associated with best solution from previous tentative level shift.

    This output is present just if input option model.lshift>0.

    Likloc

    matrix of size (2*lshiftlocref.wlength+1)-by-3 containing local sum of squares of residuals in order to decide best position of level shift:

    1st col = position of level shift;

    2nd col = local sum of squares of huberized residuals;

    3rd col = local sum of squares of raw residuals.

    This output is present just if input option model.lshift>0.

    RES

    Matrix of size T-by-(T-lshift) containing scaled residuals for all the T units of the original time series monitored in steps lshift+1, lshift+2, ..., T-lshift, where lshift+1 is the first tentative level shift position, lshift +2 is the second level shift position, and so on. This output is present just if input option model.lshift>0.

    yhat

    vector of fitted values after final (NLS=non linear least squares) step.

    $ (\hat \eta_1, \hat \eta_2, \ldots, \hat \eta_T)'$

    residuals

    Vector T-by-1 containing the scaled residuals from after final NLS step.

    weights

    Vector containing weights after adaptive reweighting. The elements of this vector are 0 or 1. These weights identify the observations which are used to compute the final NLS estimate.

    scale

    Final scale estimate of the residuals using final weights.

    \[ \hat \sigma = cor \times \sum_{i \in S_m} [y_i- \eta(x_i,\hat \beta)]^2/(m-p) \] where $S_m$ is a set of cardinality $m$ which contains the units not declared as outliers, $p$ is the total number of estimated parameters and $cor$ is a correction factor to make the estimator consistent.

    conflev

    confidence level which is used to declare outliers.

    Remark: scalar out.conflev will be used to draw the horizontal lines (confidence bands) in the plots

    outliers

    vector containing the list of the units declared as outliers using confidence level specified in input scalar conflev.

    singsub

    Number of subsets wihtout full rank. Notice that if this number is greater than 0.1*(number of subsamples) a warning is produced on the screen

    invXX

    $cov(\beta)/\hat \sigma^2$. p-by-p, square matrix.

    If the model is linear out.invXX is equal to $(X'X)^{-1}$, else out.invXX is equal to $(A'A)^{-1}$ where $A$ is the matrix of partial derivatives. More precisely:

    \[ a_{i,j}=\frac{\partial \eta_i(x_i, \hat \beta)}{\partial \hat \beta_j} \] where \begin{eqnarray} y_i & = & \eta(x_i,\beta)+ \epsilon_i \\ & = & \eta_i +\epsilon_i \\ & = & \eta(x_i,\hat \beta)+ e_i \\ & = & \hat \eta_i + e_i \end{eqnarray}

    y

    response vector y.

    X

    data matrix X containing trend, seasonal, expl and lshift, if the model is linear or linearized version of $\eta(x_i, \beta)$ if the model is non linear containing in the columns partial derivatives evaluated in correspondence of out.B(:,1) with respect to each parameter. In other words, the $i,j$-th element of out.X is \[ \frac{\partial \eta_i(x_i, \hat \beta)}{\partial \hat \beta_j} \]

    $j=1, 2, \ldots, p$, $i \in S_m$.

    The size of this matrix is:

    n-length(out.outliers)-by-p The field is present only if option yxsave is set to 1.

    class

    'LTSts'.

    varargout —Indices of the subsamples extracted for computing the estimate (the so called elemental sets) for each tentative level shift position. C : cell

    C{1} is associated with the subsamples for first tentative level shift position;

    C{2} is associated with the subsamples for second tentative level shift position;

    ...

    C{end} is associated with the subsamples for last tentative level shift position;

    References

    Rousseeuw, P.J., Perrotta D., Riani M. and Hubert, M. (2018), Robust Monitoring of Many Time Series with Application to Fraud Detection, "Econometrics and Statistics". [RPRH]

    See Also

    |

    This page has been automatically generated by our routine publishFS