LTSts extends LTS estimator to time series

It is possible to introduce a trend (up to third order), 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).

```
```

Airline data: linear trend + just one harmonic for seasonal
component.`out`

=LTSts(`y`

,
`Name, Value`

)

No seasonal component.

n=45; a=1; b=0.8; sig=1; seq=(1:n)'; y=a+b*seq+sig*randn(n,1); y(round(n/2):end)=y(round(n/2):end)+10; % model with a quadratic trend, non seasonal and level shift model=struct; model.trend=2; 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);

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')

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')

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')

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')

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')

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')

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')

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); % generate the wedgeplot % wedgeplot(out,'transpose',true,'extradata',[y1 out.yhat]);

A persisting level shift and three isolated outliers, two of which in proximity of the level shift.

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); % generate the wedgeplot % wedgeplot(out,'transpose',true,'extradata',[y1 out.yhat]);

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; % Estimate the parameters [out, varargout]=LTSts(y1,'model',model,'nsamp',500,... 'plots',2,'lshiftlocref',lshiftlocref,'msg',0); % generate the wedgeplot % wedgeplot(out,'transpose',true,'extradata',[y1 out.yhat]);

`y`

— Time series to analyze.
Vector.A row or a column vector with T elements, which contains the time series.

**
Data Types: **`single| double`

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
```

.

```
'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 interecept (default), trend = 2 implies quadratic trend ... |

`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))$. |

`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). 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; |

**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 ofrobust 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`

`out`

— description
StructureA 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 forming best initial elemental subset (that is elemental subset which produced the smallest value of the target function). |

`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(\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+lshift. 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 : cellC{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;

Rousseeuw, P.J., Perrotta D., Riani M., Hubert M. (2017), Robust Monitoring of Many Time Series with Application to Fraud Detection,

submitted.(RPRH)