Four sample programs by Roncalli:

Author: Thierry Roncalli <Thierry.Roncalli@montesquieu.u-bordeaux.fr>
Date: 2 Sep 1996
This code is provided without guarantee for public non-commercial use.


I.
Abil EL BATJI asked me how to modify the cycle2.prg
(estimation of the cycle + noise model) program of
TSM to estimate the local linear trend + cycle model.

The following program will perhaps interest other TSM users.


new;
library tsm,optmum,pgraph;

load rainfall[] = rainfall.asc;

y = rainfall[1:131];
y = y - 142.2;

/*
**  Be careful, the stationary form of the local linear
**  trend model is Delta^2(y) (and not Delta(y) which
**  corresponds to the stationary form of the random
**  walk plus noise model)
*/

y = y - lag1(y);
y = y - lag1(y);

/*
**  Using the parametrization
*/

_cycle_prmt = 0;

_tsm_parnm = "rho"|"lambda"|"sig_eta"|
             "sig_zeta"|"sig_kappa"|"sig_eps";
_tsm_optmum = 1;
_print = 1;
_tsm_Mcov = 1;
_fourier = 0;

sv = 0.98|1.6|0.1|0.01|0.01|10;

{theta,stderr,Mcov,Logl} = FD_ml(y,&sgf,sv);

/*
** After transformation, we obtain
*/

sv = theta;
sv[1] = exp(-sv[1])/(1+exp(-sv[1]));
sv[2] = pi*exp(-sv[2])/(1+exp(-sv[2]));

print "parameters  Estimated values ";
print "**********  **************** ";
print;
call printfmt(_tsm_parnm~sv,0~1);
wait;


{lambda,I} = PDGM(y);

/* Parzen lag window with bandwidth = 7 */

_smoothing = 7|5|0|0.23;
I = smoothing(I);
g = sgf(theta,lambda);
q = trunc(rows(lambda)/2);


graphset;
  _pdate = ""; _pnum = 2; fonts("simplex simgrma");
  _plegstr = "Periodogram of the Rainfall data"\
             "\000Estimated spectral generating function";
  _plegctl = {2 5 1.5 5};
  xlabel("frequency");
  xtics(0,pi,pi/4,0);
  lab = " 0 \202p\201/4 \202p\201/2 \2013\202p\201/4 \202p\201";
  asclabel(lab,0);
  xy(lambda[1:q],I[1:q,.]~g[1:q]);



/*
**  Spectral generating function for a cyclical trend model
*/

proc sgf(theta,lambda);
  local rho,lc,sig_eta,sig_zeta,sig_kappa,sig_eps;
  local theta1,theta2,g_phi,g_mu,g;

  rho = theta[1];
  lc = theta[2];
  sig_eta = theta[3];
  sig_zeta = theta[4];
  sig_kappa = theta[5];
  sig_eps = theta[6];

  /* Parameters of the spectral generating function
     of the cycle model */

  theta1 = rho|lc|sig_kappa;
  g_phi = _cycle_sgf(theta1,lambda);


  /* Parameters of the spectral generating function
     of the local linear trend model */

  theta2 = sig_eta|sig_zeta|sig_eps;
  g_mu = _cycle_sgf(theta2,lambda);


  /* Harvey's formula 2.4.60, page 64 */
  g = g_mu + 2*(1-cos(lambda)).*g_phi;

  retp(g);
endp;




II.
A user of TSM encountered some difficulties to
estimate a factor time series model, because
the estimation in time domain was sensitive to
the initialization of the Kalman filter.

A possible answer is to estimate the model in the
Frequency domain and then to estimate the factors
using the Kalman filter with diffuse prior for
initialization.

The following example may interest other TSM users.




/*

The factor time series model consists of the state space model

 (1a)  x(t) =  Z * f(t) + epsilon(t)
 (1b)  f(t) =  f(t-1) + nu(t)

with  epsilon(t)~N(0,H)   and   nu(t)~N(0,Q).


The stationary form of the state space model is

 (2)   x(t)-x(t-1) = Z*nu(t) + (1-L)*epsilon(t)

where L is the lag operator.


The corresponding spectral generating function is

 (3)   G(lambda_j) = Z*Q*Z' + 2*(1-cos(lambda_j))*H

where lambda_j is the Fourier frequency.


We can estimate the parameters of the model (1) by
maximum likelihood in the frequency domain (see
the paragraph "A note on maximum likelihood of
multivariate processes in the frequency domain"
of the TSM manual, page 32). Then, we can estimate
the factors with the Kalman filter.


In the example below, two time series x(t) = [x1(t) x2(t)]'
depend on two factors f(t) = [f1(t) f2(t)]'.

*/


new;
library tsm,optmum,pgraph;
TSMset;

/* Simulation of the model */

rndseed 123456;

Z = {1 2,3 -4};
d = {0,0};
H = {0.25 0,0 0.10};
T = eye(2);
c = {0,0};
R = eye(2);
Q = {2 0,0 1};

call SSM_build(Z,d,H,T,c,R,Q,0); /* Build the state space model */
Ns = 100;                        /* Number of simulations */
{x,f} = RND_SSM(10|-10,Ns);      /* Simulation of xt and ft */


/* Compute the spectral generating function */

proc sgf(theta,lambda);
  local Nstar,Z,Pchol,H,Q,j,gj,g;

  Nstar = rows(lambda);

  Z = reshape(theta[1:4],2,2)';
  H = (theta[5]~0)|(0~theta[6]); H = H^2;
  Q = (theta[7]~0)|(0~theta[8]); Q = Q^2;

  g = zeros(Nstar,4);

  j = 1;
  do until j > Nstar;
    gj = Z*Q*Z' + 2*(1-cos(lambda[j]))*H;
    g[j,.] = vec(gj).';
    j = j + 1;
  endo;

  retp(g);
endp;


/* Compute the multivariate periodogram */

dx = x - lag1(x);
{lambda,Ix} = PDGM2(dx);


/* Compute the log-likelihood in the frequency domain */

proc ml(theta);
  local Gx,Nstar,logl,j,Ij,Gj;

  Gx = sgf(theta,lambda);
  Nstar = rows(lambda);

  logl = zeros(Nstar,1);

  j = 1;
  do until j > Nstar;
   Ij = xpnd2(Ix,j);
   Gj = xpnd2(Gx,j);
   logl[j] = -0.5*ln(det(Gj)) -0.5*sumc(diag((inv(Gj)*Ij)));
   j = j + 1;
  endo;

  logl = real(logl);

  retp(logl);
endp;


/* Estimation of the parameters vector */


sv = vec(Z)|sqrt(diag(H)|diag(Q));
_tsm_Mcov = 1;
{theta,stderr,Mcov,Logl} = TD_ml(&ml,sv);


/* Estimation of the factors with the Kalman filter */

Z = reshape(theta[1:4],2,2)';
H = (theta[5]~0)|(0~theta[6]); H = H^2;
Q = (theta[7]~0)|(0~theta[8]); Q = Q^2;

call SSM_build(Z,d,H,T,c,R,Q,0); /* Build the state space model */
a0 = zeros(2,1);                 /* Initialize the a0 vector    */
P0 = eye(2)*1000;                /* Diffuse Prior               */
call KFiltering(x,a0,P0);        /* Run the Kalman filter       */
at = KF_matrix(3);               /* Read the state vectors a(t) */

graphset;
  _pdate = ""; _pnum = 2;
  title("True factors and estimated factors");
  _pltype = 6|6|1|1;
  xy(seqa(1,1,Ns),f~at);


parnm = 0$+"P"$+ftocv(seqa(1,1,rows(theta)),2,0);
print "                 True values    FDML estimates   ";
call printfmt(parnm~sv~theta,0~1~1);


/*********************************************************
                 True values    FDML estimates
     P01                1       0.80956994
     P02                3        3.0908726
     P03                2        1.8444238
     P04               -4       -4.0558464
     P05              0.5       0.86743126
     P06       0.31622777       0.47991068
     P07        1.4142136        1.4698517
     P08                1       0.91785868
**********************************************************/








III.
Adam Kurpiel has written a non-linearity test procedure
based on the LWG test (or the neural net test) and the
Fourier Transform surrogate.

References :

Theiler, Eubank, Longtin, Galdrikiana and Farmer [1992],
Testing for nonlinearity on time series: the method of
surrogate data, Physica D, 58, 77-94

The following program will perhaps interest other TSM users.

Th. Roncalli



new;
library tsm,optmum;
TSMset;

/* Simulate an arma(1,1) process */

{data1,retcode} = RND_arfima(0.5|-0.2|0,1,1,1,1000,512,1);

data2 = (data1 + cos(data1) )^2;

_fourier = 0;  /* Using the fast transform */
_surrogate_random = 1; /* uniform [0,2pi] randomization of the phases */

seed = 123;
call Test(data1,25,3,4,seed);

print;

seed = 123;
call Test(data2,25,3,4,seed);

/**************************************************************

----------------------------------------------
LWG-Surrogate test for neglected non-linearity
----------------------------------------------

Number of surrogate:                 25

LWG test for the original data:      0.86255
LWG test for the surrogate data:      4.31888

Surrogate test:                       0.95885
p_value:                              0.33764

----------------------------------------------
LWG-Surrogate test for neglected non-linearity
----------------------------------------------

Number of surrogate:                 25

LWG test for the original data:      10.85320
LWG test for the surrogate data:      3.52758

Surrogate test:                       2.53326
p_value:                              0.01130

***************************************************************/


/* {sigma,pvalue} = Test(y,Ns,order1,order2,seed} */

/* y: data                                      */
/* Ns: number of FT surrogate simulations       */
/* order1:  the linear autoregressive order     */
/* order2:  the non-linear autoregressive order */
/* seed: seed parameter                         */

/* sigma: the surrogate test                    */
/* pvalue: p-value of the surrogate test        */


proc (2) = Test(y,Ns,order1,order2,seed);
  local simul,TR2,pvalue,TR2s,s,ys,sigma;

  y = packr(y);

  if cols(y) /= 1;
    ERRORLOG "error: Not valid for multivariate data.";
    end;
  endif;

  {TR2,pvalue} = LWG(y,order1,order2,seed);

  simul = {};

  s = 1;
  do until s > Ns;
    cls; locate 4,5; printdos ftos(s,"surrogate N+%lf",5,0);
    ys = surrogate(y);
    simul = simul~ys;
    s = s + 1;
  endo;

  TR2s = zeros(Ns,1);

  s = 1;
  do until s > Ns;
    {TR2s[s],pvalue} = LWG(simul[.,s],order1,order2,seed);
    s = s + 1;
  endo;

  sigma = abs(TR2 - meanc(TR2s))/stdc(TR2s);
  pvalue = erfc(sigma/sqrt(2));

  if _print == 1;

    cls;
    print "----------------------------------------------";
    print "LWG-Surrogate test for neglected non-linearity";
    print "----------------------------------------------";
    print;
    print ftos(Ns,"Number of surrogate:                 %lf",1,0);
    print;
    print ftos(TR2,"LWG test for the original data:      %lf",5,5);
    print ftos(meanc(TR2s),"LWG test for the surrogate data:      %lf",5,5);
    print;
    print ftos(sigma,"Surrogate test:                       %lf",5,5);
    print ftos(pvalue,"p_value:                              %lf",5,5);

  endif;

  retp(sigma,pvalue);
endp;


proc (2) = LWG(data,p,J,seed);
  local tmin,tmax,t,Nobs,y,x,i,xi,z;
  local K,beta,u,phi,M,delta,e,R2,TR2,pvalue;

  /* target data */

  tmin = 0.1; tmax = 0.9;
  t = tmin + (tmax-tmin)*(data-minc(data)')./(maxc(data)'-minc(data)');

  Nobs = rows(t);

  /* endogeneous data */

  y = t[p+1:Nobs];

  /* exogeneous data */

  x = {};
  i = 1;
  do until i > p;
    xi = t[1+p-i:Nobs-i];
    x = x~xi;
    i = i+1;
  endo;

  /* Add a constant */

  z = ones(Nobs-p,1)~x;

  I = cols(z);    /* Input Number */
  K = cols(y);    /* Output Number */

  /* OLS: data(t) = beta0 + beta1*data(t-1) + ... + betap*data(t-p) + u(t) */

  beta = y/z;

  /* Residuals */

  u = y - z*beta;

  /* Neural Net */

  phi = Activation(z,J,seed);

  /* Form the regressor matrix */

  M = z~phi;

  /* Residuals of the Neural Net model */

  delta = u/M;
  e = u - M*delta;

  /* T*R^2 */

  R2 = 1 - (e'e)/(u'u);
  TR2 = (Nobs-p)*R2;

  /* pvalue */

  pvalue = cdfchic(TR2,J);

  retp(TR2,pvalue);
endp;


proc Net_Function(x);
  local y;
  y = 1./(1+exp(-x));
  retp(y);
endp;


proc Activation(x,J,seed);
  local Nx,u,ht,at;

  rndseed seed;

  Nx = cols(x);

  u = 4*rndu(J,Nx) - 2;  /* [-2,2] uniform random */
  ht = x*u';
  at = Net_Function(ht);

  retp(at);
endp;


Iv.
Common/specific component model

In the paper "The second bank of the United States in
its early years: lessons from a transition towards an
inter-regional monetary union.", M-A. Senegas is
interested in the convergence process of five american
regional inflation rates at the beginning of the
nineteenth century (1817 - 1840). A way to assess the
evolution of the convergence degree in these rates goes
through the estimation of their potential common trend
throughout the whole period this aim, a common/specific
component model is applied and estimated by Kalman filtering
in the time domain. Each regional inflation rate is
assumed to be determined, on the one hand by a component
which is common to all regional inflation rates
in the Union and, on the other hand, by an indiosyncratic
shock specific to each of them. The two (common and
specific) components are supposed to be I(1) processes

M-A. Senegas uses the following state-space form

  (1a)   y(t) = Z * a(t)
  (1b)   a(t) = a(t-1) + nu(t)
where
         nu(t)~N(0,Q)
and
         Q = diag(sig_c^2,sig1^2,sig2^2,sig3^2,sig4^2,sig5^2)

y(t) = (pi1(t),pi2(t),...,pi5(t))' where pii(t) refers to the
inflation rate of region i (New-York, Philadelphia, Charleston,
New Orleans and Cincinnati, source: Cole, A. (1938), Wholesale
commodity prices in the United States, 1700-1861, Harvard U. Press).
a(t) = (C(t),mu1(t),...,mu5(t))' where C(t) is the common component
and mui(t) is the specific component relative to region i.

The measurement equation has the following form:
   _        _    _                 _   _        _
   | pi1(t) |    | beta1 1 0 0 0 0 |   | C(t)   |
   | pi2(t) |    | beta2 0 1 0 0 0 |   | mu1(t) |
   | pi3(t) | =  | beta3 0 0 1 0 0 | * | mu2(t) |
   | pi4(t) |    | beta4 0 0 0 1 0 |   | mu3(t) |
   | pi5(t) |    | beta5 0 0 0 0 1 |   | mu4(t) |
   -        -    -                 -   | mu5(t) |
                                       -        -
M-A Senegas has estimated his model in time domain. Note that
we can also estimate the parameters in the frequency domain.
The stationary form is y(t)-y(t-1) = Z*nu(t). The corresponding
spectral generating function is G(lambda_j) = Z*Q*Z'.

The following programs will perhaps give some ideas
to other TSM users. This article is available upon request
to the author. Email:senegas@montesquieu.u-bordeaux.fr


new;
library tsm,optmum,pgraph;
TSMset;

load y = pigen;

y = packr(y);

d = zeros(5,1);
H = zeros(5,5);
T = eye(6);
c = zeros(6,1);
R = eye(6);

a0 = 0|y[1,.]';
P0 = zeros(6,6);

/* Compute the log-likelihood in the time domain */

proc ml(theta);
  local beta,sigma,Z,Q,logl;

  beta = theta[1:5];
  sigma = theta[6:11]^2;

  Z = beta~eye(5);
  Q = diagrv(zeros(6,6),sigma);

  call SSM_build(Z,d,H,T,c,R,Q,0); /* Build the state space model */
  call KFiltering(y,a0,P0);        /* Run the Kalman filter       */
  logl = KF_ml;              /* Compute the log-likelihood vector */

  retp(logl);
endp;


/* Estimation of the parameters vector */

sv = ones(11,1);
_tsm_Mcov = 1;
{theta,stderr,Mcov,Logl} = TD_ml(&ml,sv);

load theta;

beta = theta[1:5];
sigma = theta[6:11]^2;

Z = beta~eye(5);
Q = diagrv(zeros(6,6),sigma);


call SSM_build(Z,d,H,T,c,R,Q,0); /* Build the state space model */
call KFiltering(y,a0,P0);        /* Run the Kalman filter       */
at = KF_matrix(3);               /* Read the state vectors a(t) */

/* Read the covariance matrices of a(t) */

Pt = KF_matrix(4);

/* Common component +/- 1.96 * SE */

bound = at[.,1] + (1.96~-1.96).*sqrt(Pt[.,1]);

graphset;
  _pdate = ""; _pnum = 2;
  lab=" [1817] [1818] [1819] [1820] [1821] "\
      " [1822] [1823] [1824] [1825] [1826] [1827] [1828] "\
      " [1829] [1830] [1831] [1832] [1833] [1834] [1835] "\
      " [1836] [1837] [1838] [1839] [1840] [1841] ";
  asclabel(lab,0);
  xtics(1,300,12,4);
  title("Common component significance"\
        "\Lregional inflation rates"\
        "\L(1817-1841)");

  fonts("simplex,simgrma");

  _plegstr="C]t|t-1[\000C]t|t-1[-1.96*SE\000C]t|t-1[+1.96*SE";
  _plegctl = {2 5 1 0.7};

  _pltype = {6 5 5};
  _plwidth = 5|5|5;
  xy(seqa(1,1,rows(at)),at[.,1]~bound);


/*
** Here is the same program, but the estimation is performed
** in the Frequency domain
*/

/*
new;
library tsm,optmum,pgraph;
TSMset;

load y = pigen;

y = packr(y);

dy = y - lag1(y);

_fourier = 1;
{lambda,Iy} = PDGM2(dy);

/* Compute the spectral generating function */

proc sgf(theta,lambda);
  local Nstar,beta,sigma,Z,Q,j,gj,g;

  Nstar = rows(lambda);

  beta = theta[1:5];
  sigma = theta[6:11]^2;

  Q = diagrv(zeros(6,6),sigma);
  Z = beta~eye(5);

  g = vec(Z*Q*Z').' + zeros(Nstar,1);

  retp(g);
endp;


/* Compute the log-likelihood in the frequency domain */

proc ml(theta);
  local Gy,Nstar,logl,j,Ij,Gj;

  Gy = sgf(theta,lambda);
  Nstar = rows(lambda);

  logl = zeros(Nstar,1);

  j = 1;
  do until j > Nstar;
   Ij = xpnd2(Iy,j);
   Gj = xpnd2(Gy,j);
   logl[j] = -0.5*ln(det(Gj)) -0.5*sumc(diag((inv(Gj)*Ij)));
   j = j + 1;
  endo;

  logl = real(logl);

  retp(logl);
endp;


/* Estimation of the parameters vector */


sv = ones(11,1);
_tsm_Mcov = 1;
{theta,stderr,Mcov,Logl} = TD_ml(&ml,sv);

*/


