From pgilbert at bank-banque-canada.ca  Fri Sep 14 17:50:35 2007
From: pgilbert at bank-banque-canada.ca (Paul Gilbert)
Date: Fri, 14 Sep 2007 11:50:35 -0400
Subject: [Rd] partial correlation function for multivariate time series
In-Reply-To: <3c12769c0709140748o137cb1ebg66f33d9790014c4d@mail.gmail.com>
References: <3c12769c0709100805r4ca5785bq3a932665c82084fd@mail.gmail.com>	
	<46E6BBCC.5040907@bank-banque-canada.ca>
	<3c12769c0709140748o137cb1ebg66f33d9790014c4d@mail.gmail.com>
Message-ID: <46EAADCB.20305@bank-banque-canada.ca>

Simone Giannerini wrote:
> Dear Paul,
> 
> there is no mention to the pacf in the multivariate setting in the book 
> you suggested.

My apologies, I should have looked more carefully and realized the pacf 
discussion in Granger and Newbold is all univariate.

> As I mentioned in private I suspect that pacf() in the multivariate case 
> computes the
> partial autoregression matrix (in the terminology of Reinsel) rather 
> than the partial autocorrelation matrix

It looks like pacf() uses the calculation in stats:::ar.yw.default , so 
roughly sounds consistent with what you say. I do think parts of this 
code pre-date Reisel's book, so inconsistency with his book would 
probably be imperfect foresight.

BTW, I think bug reports with suggested fixes are usually more useful. 
And, at the very least, it seems the references in the documentation 
could be improved.

Paul Gilbert

> as the two coincide in the univariate case but are different in the 
> multivariate case as stated explicitly
> in Reisel (Sec. 3.3).
> This would explain the coefficients well above 1 (in modulus) in the 
> example I have given.
> To support my statement you can fit A VAR model to the data and compare 
> the coefficients with the results
> from pacf():
> 
> set.seed(10)
> x <- rnorm(1000,sd=10000)
> y <- rnorm(1000,sd=1)
> 
> library(vars);
> mod1 <- VAR(ts(cbind(x,y)),p=4,type="none"); # fit a VAR (OLS)
> 
> ## Have a look at estimated coefficients
> 
>>  noquote(formatC(mod1$varresult$x$coefficients,format="f"))
> ## Edited by me, compare with the first column of the results below from 
> the pacf (pacf(....)$acf[,1,])
>      x.l1     x.l2    x.l3    x.l4   
>      0.047    0.013   0.004   0.014
>      y.l1     y.l2    y.l3    y.l4  
>      374.117  72.758 -526.452 126.610
>>  noquote(formatC(mod1$varresult$y$coefficients,format="f"))
> ## Edited by me, compare with the second column of the results below 
> from the pacf (pacf(....)$acf[,2,])
>      x.l1    x.l2     x.l3    x.l4 
>      0.000   -0.000   0.000   0.000 
>      y.l1    y.l2     y.l3    y.l4
>     -0.046   -0.025  -0.033  -0.020
> 
> pacf(ts(cbind(x,y)),plot=FALSE,lag.max=4)
> 
> Partial autocorrelations of series 'ts(cbind(x, y))', by lag
> 
> , , x
> 
>  x             y           
>     0.047 ( 1)    0.000 (-1)
>     0.011 ( 2)    0.000 (-2)
>     0.005 ( 3)    0.000 (-3)
>     0.013 ( 4)    0.000 (-4)
> 
> , , y
> 
>  x             y           
>   374.052 ( 1)   -0.045 ( 1)
>    66.717 ( 2)   -0.024 ( 2)
>  -535.810 ( 3)   -0.031 ( 3)
>   120.802 ( 4)   -0.020 ( 4)
>  
> As you can see the coefficients fairly agree.
> I might file a bug report in some days unless someone will prove me 
> wrong before.
> 
> Regards,
> 
> Simone
> 
> 
> 
> On 11/09/2007, *Paul Gilbert* <pgilbert at bank-banque-canada.ca 
> <mailto:pgilbert at bank-banque-canada.ca>> wrote:
> 
>     I think the reference for pacf is
> 
>     @BOOK{GraNew77,
>        author =    {Granger, C. W. J. and Newbold, Paul},
>        title =     {Forecasting Economic Time Series},
>        publisher = {Academic Press},
>        year =      1977
>        }
> 
>     It certainly would not be Reisel's book, as parts of the code predate
>     that by many years.
> 
>     Paul Gilbert
> 
>     Simone Giannerini wrote:
>      > Dear all,
>      >
>      > I found the following behaviour with pacf() in the multivariate
>     case,
>      >
>      > set.seed(10)
>      > x <- rnorm(1000,sd=10000)
>      > y <- rnorm(1000,sd=1)
>      > pacf(ts(cbind(x,y)),plot=FALSE,lag.max=10)
>      >
>      > Partial autocorrelations of series 'cbind(x, y)', by lag
>      >
>      > , , x
>      >
>      >  x              y
>      >     0.047 (  1)    0.000 ( -1)
>      >     0.011 (  2)    0.000 ( -2)
>      >     0.005 (  3)    0.000 ( -3)
>      >     0.013 (  4)    0.000 ( -4)
>      >     0.050 (  5)    0.000 ( -5)
>      >     0.034 (  6)    0.000 ( -6)
>      >     0.026 (  7)    0.000 ( -7)
>      >    -0.029 (  8)    0.000 ( -8)
>      >    -0.010 (  9)    0.000 ( -9)
>      >    -0.013 ( 10)    0.000 (-10)
>      >
>      > , , y
>      >
>      >  x              y
>      >   374.052 (  1)   -0.045 (  1)
>      >    66.717 (  2)   -0.024 (  2)
>      >  -535.810 (  3)   -0.031 (  3)
>      >   120.802 (  4)   -0.020 (  4)
>      >   142.824 (  5)    0.004 (  5)
>      >  -211.711 (  6)   -0.010 (  6)
>      >   201.856 (  7)    0.058 (  7)
>      >   286.079 (  8)   -0.035 (  8)
>      >  -134.057 (  9)    0.032 (  9)
>      >   -18.200 ( 10)    0.019 ( 10)
>      >
>      > Since there are multiple ways of defining the pacf for
>     multivariate time
>      > series
>      > (see e.g. G.C. Reinsel, Elements of multivariate time series
>     analysis, II
>      > edition, Springer, section 3.3) and given that
>      > in ?pacf there is no reference to articles or books regarding its
>      > computation
>      > I do not know whether this behaviour is expected or it is a bug
>     instead.
>      > In the first case could you provide a reference for it? In the
>     second case I
>      > might file a bug report.
>      > Thank you for the great work you are doing for the scientific
>     community.
>      >
>      > kind regards,
>      >
>      > Simone Giannerini
>      >
>      > WINDOWS
>      >
>      > platform       i386-pc-mingw32
>      > arch           i386
>      > os             mingw32
>      > system         i386, mingw32
>      > status
>      > major          2
>      > minor          5.1
>      > year           2007
>      > month          06
>      > day            27
>      > svn rev        42083
>      > language       R
>      > version.string R version 2.5.1 (2007-06-27)
>      >
>      > LINUX
>      >
>      >
>      >>R.Version()
>      >
>      > $platform
>      > [1] "x86_64-unknown-linux-gnu"
>      >
>      > $arch
>      > [1] "x86_64"
>      >
>      > $os
>      > [1] "linux-gnu"
>      >
>      > $system
>      > [1] "x86_64, linux-gnu"
>      >
>      > $status
>      > [1] ""
>      >
>      > $major
>      > [1] "2"
>      >
>      > $minor
>      > [1] "5.1"
>      >
>      > $year
>      > [1] "2007"
>      >
>      > $month
>      > [1] "06"
>      >
>      > $day
>      > [1] "27"
>      >
>      > $`svn rev`
>      > [1] "42083"
>      >
>      > $language
>      > [1] "R"
>      >
>      > $version.string
>      > [1] "R version 2.5.1 (2007-06-27)"
>     ====================================================================================
> 
>     La version fran?aise suit le texte anglais.
> 
>     ------------------------------------------------------------------------------------
> 
> 
>     This email may contain privileged and/or confidential information,
>     and the Bank of
>     Canada does not waive any related rights. Any distribution, use, or
>     copying of this
>     email or the information it contains by other than the intended
>     recipient is
>     unauthorized. If you received this email in error please delete it
>     immediately from
>     your system and notify the sender promptly by email that you have
>     done so.
> 
>     ------------------------------------------------------------------------------------
> 
> 
>     Le pr?sent courriel peut contenir de l'information privil?gi?e ou
>     confidentielle.
>     La Banque du Canada ne renonce pas aux droits qui s'y rapportent.
>     Toute diffusion,
>     utilisation ou copie de ce courriel ou des renseignements qu'il
>     contient par une
>     personne autre que le ou les destinataires d?sign?s est interdite.
>     Si vous recevez
>     ce courriel par erreur, veuillez le supprimer imm?diatement et
>     envoyer sans d?lai ?
>     l'exp?diteur un message ?lectronique pour l'aviser que vous avez
>     ?limin? de votre
>     ordinateur toute copie du courriel re?u.
> 
> 
> 
> 
> -- 
> ______________________________________________________
> 
> Simone Giannerini
> Dipartimento di Scienze Statistiche "Paolo Fortunati"
> Universita' di Bologna
> Via delle belle arti 41 - 40126  Bologna,  ITALY
> Tel: +39 051 2098262  Fax: +39 051 232153
> http://www2.stat.unibo.it/giannerini/
> ______________________________________________________
====================================================================================

La version fran?aise suit le texte anglais.

------------------------------------------------------------------------------------

This email may contain privileged and/or confidential inform...{{dropped}}



[Rd] partial correlation function for multivariate time series
Paul Gilbert