[R] Comparing complex numbers

2008-07-11 Thread David Stoffer

Is there an easy way to compare complex numbers?

Here is a small example: 

>   (z1=polyroot(c(1,-.4,-.45)))
[1]  1.11-0i -2.00+0i
>   (z2=polyroot(c(1,1,.25)))
[1] -2+0i -2+0i
>   x=0
>   if(any(identical(z1,z2))) x=99
>x
[1] 0

#  real and imaginary parts:

>Re(z1); Im(z1)
[1]  1.11 -2.00
[1] -8.4968e-21   8.4968e-21

>Re(z2); Im(z2)
[1] -2 -2
[1]  0  0

Both z1 and z2 have a root of -2, but I guess Im(z1)
isn't close enough to zero for identical().

TIA.


-
The power of accurate observation is commonly called cynicism 
by those who have not got it.  George Bernard Shaw
-- 
View this message in context: 
http://www.nabble.com/Comparing-complex-numbers-tp18406924p18406924.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Comparing complex numbers

2008-07-11 Thread David Stoffer

Thanks- I was hoping I wouldn't have to loop, but using abs() is better than 
comparing Re()s and Im()s, which was my original thought.  I have to compare
all the roots, so I used something like this:

> for (i in 1:length(z1))  if(any(abs(z1[i]-z2[1:length(z2)])<1e-15))
> print("ouch")

... thanks again for your help.



Duncan Murdoch-2 wrote:
> 
> On 7/11/2008 11:51 AM, David Stoffer wrote:
>> Is there an easy way to compare complex numbers?
>> 
>> Here is a small example: 
>> 
>>>   (z1=polyroot(c(1,-.4,-.45)))
>> [1]  1.11-0i -2.00+0i
>>>   (z2=polyroot(c(1,1,.25)))
>> [1] -2+0i -2+0i
>>>   x=0
>>>   if(any(identical(z1,z2))) x=99
>>>x
>> [1] 0
>> 
>> #  real and imaginary parts:
>> 
>>>Re(z1); Im(z1)
>> [1]  1.11 -2.00
>> [1] -8.4968e-21   8.4968e-21
>> 
>>>Re(z2); Im(z2)
>> [1] -2 -2
>> [1]  0  0
>> 
>> Both z1 and z2 have a root of -2, but I guess Im(z1)
>> isn't close enough to zero for identical().
> 
> == is the test you had in mind, I think, but like identical it requires 
> exact equality.  identical() wants the whole vector to be identical. 
> all.equal() checks for element by element approximate equality.
> 
> So you can get element by element approximate equality by something like
> 
>  > for (i in 1:length(z1)) print(all.equal(z1[i], z2[i]))
> [1] "Mean relative Mod difference: 2.8"
> [1] TRUE
> 
> The trouble with this approach is that it will use a different scale for 
> each comparison.  I think you really need to fashion your own test, e.g.
> 
>  > abs(z1-z2) < 1e-15
> [1] FALSE  TRUE
> 
> Duncan Murdoch
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 


-
The power of accurate observation is commonly called cynicism 
by those who have not got it.  George Bernard Shaw
-- 
View this message in context: 
http://www.nabble.com/Comparing-complex-numbers-tp18406924p18411302.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Box.test degrees of freedom

2008-08-09 Thread David Stoffer

I believe tsdiag() uses the correct degrees of freedom in applying Box.test,
but the graphic shows "lag" on the horizontal axis when it should display
"degrees of freedom".   



raf.rossignol wrote:
> 
> Hello,
> 
> Prof Brian Ripley wrote:
>> 
>> I think you are referring to its application to the residuals of an 
>> ARMA(p, q) fit, and that is not what Box.test says it does.
>> 
>> It is very easy to edit the code if you want to use a different degrees
>> of
>> freedom.
>> 
> I am also new to R, but it seems to me that there is still something
> confusing, not in Box.test but in tsdiag.Arima
> Indeed, the help of tsdiag says "The methods for 'arima' and 'StructTS'
> [...] use the Ljung-Box version of the portmanteau test." 
> So we could expect the degrees of freedom 'h-p-q' to be used, but a look
> at tsdiag.Arima shows it uses Box.test at lags h=1:gof.lag, with degrees
> of freedom equal to h, and  not h-p-q. Do you think this is a mistake in
> tsdiag.Arima or is there some experimental (or theoretical) reason
> supporting this choice ?
> Best regards
> 
> 
> 


-
The power of accurate observation is commonly called cynicism 
by those who have not got it.  George Bernard Shaw
-- 
View this message in context: 
http://www.nabble.com/Box.test-degrees-of-freedom-tp17277964p18907921.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Box.test degrees of freedom

2008-08-10 Thread David Stoffer

I stand corrected.  I thought I checked this a long time ago, but apparently
not.  tsdiag.Arima DOES NOT use the fact that the series it is testing (or
diagnosing, if you will) are residuals from an ARIMA fit. 

I keep a list of R time series bloopers here: 
http://www.stat.pitt.edu/stoffer/tsa2/Rissues.htm along with some
work-arounds over here: http://www.stat.pitt.edu/stoffer/tsa2/Examples.htm




David Stoffer wrote:
> 
> I believe tsdiag() uses the correct degrees of freedom in applying
> Box.test, but the graphic shows "lag" on the horizontal axis when it
> should display "degrees of freedom".   
> 
> 
> 
> raf.rossignol wrote:
>> 
>> Hello,
>> 
>> Prof Brian Ripley wrote:
>>> 
>>> I think you are referring to its application to the residuals of an 
>>> ARMA(p, q) fit, and that is not what Box.test says it does.
>>> 
>>> It is very easy to edit the code if you want to use a different degrees
>>> of
>>> freedom.
>>> 
>> I am also new to R, but it seems to me that there is still something
>> confusing, not in Box.test but in tsdiag.Arima
>> Indeed, the help of tsdiag says "The methods for 'arima' and 'StructTS'
>> [...] use the Ljung-Box version of the portmanteau test." 
>> So we could expect the degrees of freedom 'h-p-q' to be used, but a look
>> at tsdiag.Arima shows it uses Box.test at lags h=1:gof.lag, with degrees
>> of freedom equal to h, and  not h-p-q. Do you think this is a mistake in
>> tsdiag.Arima or is there some experimental (or theoretical) reason
>> supporting this choice ?
>> Best regards
>> 
>> 
>> 
> 
> 


-
The power of accurate observation is commonly called cynicism 
by those who have not got it.  George Bernard Shaw
-- 
View this message in context: 
http://www.nabble.com/Box.test-degrees-of-freedom-tp17277964p18917747.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Box.test degrees of freedom

2008-08-11 Thread David Stoffer

No, the df are H-(p+q), against what your intuition tells you.  The problem
is that you're thinking of the residuals as losing a df for each parameter,
but the asymptotics of the problem involve the sample autocorrelations of
the residuals.  You can read the details in the original paper by Box and
Pierce (?Box.test for the reference).   By the way, you're not alone,
Minitab makes the same mistake you did.



raf.rossignol wrote:
> 
> 
> 
> David Stoffer wrote:
>> 
>> I stand corrected.  I thought I checked this a long time ago, but
>> apparently not.  tsdiag.Arima DOES NOT use the fact that the series it is
>> testing (or diagnosing, if you will) are residuals from an ARIMA fit. 
>> 
>> I keep a list of R time series bloopers here: 
>> http://www.stat.pitt.edu/stoffer/tsa2/Rissues.htm along with some
>> work-arounds over here:
>> http://www.stat.pitt.edu/stoffer/tsa2/Examples.htm
>> 
>> 
> Thanks,
> 
> by the way, if an intercept is included in the model (which is the default
> setting for arima()), it seems to me that the good number of degrees of
> freedom sould be (h-p-q-1). Do you agree ?
> 
> 


-
The power of accurate observation is commonly called cynicism 
by those who have not got it.  George Bernard Shaw
-- 
View this message in context: 
http://www.nabble.com/Box.test-degrees-of-freedom-tp17277964p18926678.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] arima and xreg

2008-09-10 Thread David Stoffer

Your model is close, but not correct... there are no t's on the parameters
and the U's aren't lagged.  

You can find an ARMAX example on our "quick fix" page:
http://www.stat.pitt.edu/stoffer/tsa2/R_time_series_quick_fix.htm .  The
example is near the bottom and just above the spectral analysis example, but
you may want to read the "regression with autocorrelated errors" first to
get some background.  




Jose Capco wrote:
> 
> Dear R-help-archive..
> 
> I am trying to figure out how to make arima prediction when I have a
> process involving multivariate time series input, and one output time
> series (output is to be predicted) .. (thus strictly speaking its an
> ARMAX process).  I know that the arima function of R was not designed
> to handle multivariate analysis (there is dse but it doesnt handle
> arma multivariate analysis, only simulations). But there is this
> beautiful "xreg" as parameter for arima and I was wondering..
> for the case of one output series I can actually "trick" R in doing
> multivariate time series for me no?.. because I saw in the
> documentation, xreg can be inputed as a ---matrix--- with output.len
> (length of output data) number of rows.. So in fact I can let the
> different columns of xreg to actually be the different input time
> series I need!
> 
> Is anyone familiar in how arima with xreg as given estimate models? ..
> how is the model assumed?
> 
> supposing I write :
> 
> arima(y, xreg=U, order=c(3,0,2))
> 
> how is y_t calculated? (supposing U has 2 columns, with U[1] being
> first column and U[2] second column)
> 
> is it
> 
> y_t = theta_(t-1)y_t-1 +  + theta_t-3 y_t-3 + intercept + U[1]_t +
> psi[1]_t-1 U[1]_t-1 + psi[1]_t-2 U[1]_t-2 + +  psi[2]U[2]_t-2 +
> e_t + phi_t-1 e_t-1 + phi_t-2 e_t-2
> 
> ??
> 
> e_t .. etc. are the white noise series of the model.
> 
> the documentation is totally vague when it comes to xreg. I hope it is
> like above :)
> 
> Would appreciate any remarks or comments. Thanks in advance.
> 
> Sincerely,
> Jose
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 


-
The power of accurate observation is commonly called cynicism 
by those who have not got it.  George Bernard Shaw
-- 
View this message in context: 
http://www.nabble.com/arima-and-xreg-tp19415386p19427576.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] arima and xreg

2008-09-11 Thread David Stoffer

You can have lagged inputs in the xreg statement, you just have to construct
the input matrix properly so the dimensions are the same, e.g.,

x = ts.intersect(mort, trend, part, lag(part,-4))
arima(x[,1],order=c(2,0,1), xreg=x[,2:4])

... and yes you have to worry about singularities or even multicolinearity
(near or computational singularity), e.g., this fails:

x = ts.intersect(mort, trend, part, part)
arima(x[,1],order=c(2,0,1), xreg=x[,2:4])
 



Jose Capco wrote:
> 
> 
> 
> On Sep 11, 6:24 am, David Stoffer <[EMAIL PROTECTED]> wrote:
>> Your model is close, but not correct... there are no t's on the
>> parameters
>> and the U's aren't lagged.
>>
>> You can find an ARMAX example on our "quick fix"
>> page:http://www.stat.pitt.edu/stoffer/tsa2/R_time_series_quick_fix.htm. 
>> The
>> example is near the bottom and just above the spectral analysis example,
>> but
>> you may want to read the "regression with autocorrelated errors" example
>> first to get some background.
>>
>>
>>
> 
> Ok. so arima of R can only deal with unlagged inputs (thus xreg has
> the latest value in the equation). Your example was an ARX (no moving
> averages)
> here is the code of yoru example : armax.fit = arima(mort,
> order=c(2,0,0), xreg=cbind(trend, part))
> I guess I can change it to ARMAX if I use order=c(p,0,q)  =)
> 
> Now theres one thing that might be worth mentioning here though. The
> above can only work (I am guessing) if the input values (the columns
> of xreg) are uncorrelated (or what word do use for that, sorry Im a
> pure mathematician, not a statistician :p ). What I mean is that the
> matrix in which optim of R must be singular (otherwise I think, from
> my last try when using two equal valued columns in xreg, arima will
> complain that optim returns an infinity value). Is there a way to
> check if the xreg matrix have uncorrelated inputs and then just
> discard the column until xreg becomes uncorrelated? I'll do a few more
> experiment with xreg and report here.. Im doing this because the
> documentation (as many of us already know.. ) do not really explain
> xreg very well.
> 
> Sincerely,
> Jose Capco
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 


-
The power of accurate observation is commonly called cynicism 
by those who have not got it.  George Bernard Shaw
-- 
View this message in context: 
http://www.nabble.com/arima-and-xreg-tp19415386p19440600.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Help with 'spectrum'

2008-09-12 Thread David Stoffer

Kevin- this is a simple rescaling of the axes so that the "area under the
curve" remains constant (and is half of the variance since you only look at
the positive frequencies).   In this case, freq(x) = 1/dx, where dx is the
time between points.   It is basically a graphic device so that you get
pretty graphics and it's akin to drawing probability histograms so that area
corresponds to probability.  I think you'd get a good idea of what is going
on by doing this:

x <- ts(cos(2*pi*1:60*1/12))# think of monthly high temps
y <- ts(x, freq=12)
par(mfrow=c(2,1))
plot(x)
plot(y)
# see the difference?  ... which do you prefer? now do this:

par(mfrow=c(2,1))
spec.pgram(x, taper=0, log="no")# freq axis is 0 to .5 cycles per unit
time (which is 1)
spec.pgram(y, taper=0, log="no")# freq axis is 0 to 6 cycles per unit
time (which is 12)
#   ... which do you prefer?

and you'll see by stretching the frequency axis you have to adjust the
spectrum axis accordingly so that you maintain "variance under the curve".





rkevinburton wrote:
> 
> For the command 'spectrum' I read:
> 
> The spectrum here is defined with scaling 1/frequency(x), following
> S-PLUS. This makes the spectral density a density over the range
> (-frequency(x)/2, +frequency(x)/2], whereas a more common scaling is 2π
> and range (-0.5, 0.5] (e.g., Bloomfield) or 1 and range (-π, π]. 
> 
> 
> Forgive my ignorance but I am having a hard time interpreting this. Does
> this mean that in the spectrum output every element of the $spec array is
> scaled by 1/frequency(x)? I am having a hard time determing what is meant
> by 'frequency'.Say I define a time series for a year with samples for
> every day. I input a 'frequency' of 365 (which in my mind is the period).
> On the output of 'spectrum' would this mean that every element of the
> $spec array is scaled by 1/365? There is a corresponding frequency array
> on the output from 'spectrum'. If the frequency is 365 and an element in
> the frequency array output from 'spectrum' is .1 am I to assume that the
> period is 36.5 and a corresponding sin wave would be sin(2 * pi *
> 36.5/365)?
> 
> Thank you in advance for helping me clear up some confusion.
> 
> Kevin
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 


-
The power of accurate observation is commonly called cynicism 
by those who have not got it.  George Bernard Shaw
-- 
View this message in context: 
http://www.nabble.com/Help-with-%27spectrum%27-tp19396471p19466898.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] cross-correlation lag.plot?

2008-05-07 Thread David Stoffer

http://www.stat.pitt.edu/stoffer/tsa2/Examples.htm


tom soyer wrote:
> 
> Hi,
> 
> Does anyone know if R has a function that is similar to lag.plot but
> instead
> of auto-correlation, it plots cross-correlation with lags?
> 
> Thanks,
> 
> -- 
> Tom
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 


-
The power of accurate observation is commonly called cynicism 
by those who have not got it.  George Bernard Shaw
-- 
View this message in context: 
http://www.nabble.com/cross-correlation-lag.plot--tp17114046p17117624.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] estimate phase shift between two signals

2008-06-05 Thread David Stoffer

help(spec.pgram) - then look at the examples at the bottom of the page



Dylan Beaudette-3 wrote:
> 
> Hi,
> 
> Are there any functions in R that could be used to estimate the
> phase-shift 
> between two semi-sinusoidal vectors? Here is what I have tried so far,
> using 
> the spectrum() function -- possibly incorrectly:
> 
> 
> # generate some fake data, normalized to unit circle
> x <- jitter(seq(-2*pi, 2*pi, by=0.1), amount=pi/8)
> 
> # functions defining two out-of-phase phenomena
> f1 <- function(x) jitter(sin(x), amount=0.25)
> f2 <- function(x, a) jitter(sin(x + a), amount=0.25)
> 
> # compute y-values
> # we are setting the phase shift arbitrarily
> s <- pi/1.5632198
> y1 <- f1(x)
> y2 <- f2(x, s)
> 
> 
> # plot:
> plot(x, y1, type='p', col='red', cex=0.5)
> lines(lowess(x, y1, f=0.25), col='red')
> 
> points(x, y2, col='blue', cex=0.5)
> lines(lowess(x, y2, f=0.25), col='blue')
> 
> 
> # generate time series object
> comb.ts <- ts(matrix(c(y1, y2), ncol=2))
> 
> # multivariate spectral decomposition
> spec <- spectrum(comb.ts, detrend=FALSE)
> 
> # but how to interpret the phase estimate?
> mean(spec$phase)
> 
> the mean 'phase' as returned from spectrum() does not seem to match the
> value 
> used to generate the data... Am I mis-interpreting the use or output from 
> spectrum() here? If so, is there a general procedure for estimating a 
> phase-shift between two noisy signals? Would I first have to fit a smooth 
> function in order to solve this analytically?
> 
> Thanks in advance,
> 
> 
> 
> -- 
> Dylan Beaudette
> Soil Resource Laboratory
> http://casoilresource.lawr.ucdavis.edu/
> University of California at Davis
> 530.754.7341
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 


-
The power of accurate observation is commonly called cynicism 
by those who have not got it.  George Bernard Shaw
-- 
View this message in context: 
http://www.nabble.com/estimate-phase-shift-between-two-signals-tp17653636p17682957.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Kalman Filter

2008-03-04 Thread David Stoffer

Vladimir- there are at least 3 packages that will facilitate state space
modeling: 
http://cran.r-project.org/src/contrib/Descriptions/dlm.html DLM ,
http://cran.r-project.org/src/contrib/Descriptions/dse.html DSE , and
http://cran.r-project.org/src/contrib/Descriptions/sspir.html SSPIR .

In addition, I have scripts and examples for fitting state space models and
running
the Kalman filter and smoother for Chapter 6 of our text,  
http://www.stat.pitt.edu/stoffer/tsa2/ tsa2 .  Go to "R CODE (Ch 6)" using
the
blue bar at the top.  There you will find the scripts and examples.  If you
have any
questions, feel free to contact me.


Vladimír Šamaj wrote:
> 
> Hi
> 
> My name is Vladimir Samaj. I am a student of Univerzity of Zilina. I am
> trying to implement Kalman Filter into my school work. I have some
> problems
> with understanding of R version of Kalman Filter in package stats(
> functions
> KalmanLike, KalmanRun, KalmanSmooth,KalmanForecast).
> 
> 1) Can you tell me how are you seting the initial values of state vector
> in
> Kalman Filter? Are you using some method?
> 
> 2) I have fond function StructTS in stats package. I dont understand, how
> exactly, are you computing(what method are you using) fitted values which
> are the output of this function( $fitted ) . In description od this
> function
> is that it fit a structural model for a time series by maximum likehood.
> Does it means, that the fitted values are fit by maximum likehood? If so
> how
> does look the likehood function?
> 
> 3)Finaly, I dont understand smooting problem.  What I know is that, if I
> have t observations of some time serie, I can use  function KalmanRun to
> get
> estimates of state vector. And if I gain aditional observations of time
> serie( T > t ), I shoud use KalmanSmooth function to smooth estimates of
> state vector. I dont understand, that how shoud I "tell" to KalmanSmooth
> funtion that I allready did filtering and it shoud use the values from
> filtering to smoothing.
> 
> I will be glad if you help me. I hope that my folmulations were correct.
> 
> Thank you very much.
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 


-
The power of accurate observation is commonly called cynicism 
by those who have not got it.  George Bernard Shaw
-- 
View this message in context: 
http://www.nabble.com/Kalman-Filter-tp15696135p15826374.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Need help for calculating cross-correlation between 4 multivariate time series data

2008-03-05 Thread David Stoffer

You can use acf(), but it will be messy and the labeling of the plots
is confusing and perhaps misleading... check out issue number 4 
at http://www.stat.pitt.edu/stoffer/tsa2/Rissues.htm

I would recommend setting up a grid of the ccfs and you could 
automate this (i.e., write a loop)  if need be.



sandeep jose joseph wrote:
> 
> Hi all,
>  I would like to know whether there is any function in R were i
> can
> find the cross-correlation of two or more multivariate (time series) data.
> I
> tried the function ccf() but it seems like to have two univariate
> datasets.
> Please let me know.
> sincerely,
> sandeep
> 
> -- 
> Sandeep Joseph PhD
> Post Doctoral Associate
> Center for Tropical & Emerging Global Diseases
> Paul D. Coverdell Center, Rm 335A
> 500 D. W. Brooks Drive
> Athens, GA 30602-7394
> ph 404-578-6790
> 
> 2091 south Miledge Ave
> Apt A6, Athens,
> GA-30605.
> phone: 706-850-0056
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 


-
The power of accurate observation is commonly called cynicism 
by those who have not got it.  George Bernard Shaw
-- 
View this message in context: 
http://www.nabble.com/Need-help-for-calculating-cross-correlation-between-4-multivariate-time-series-data-tp15853171p15857688.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] array in version 2.8.0

2008-11-02 Thread David Stoffer

What happened?  TIA.

In version 2.7.x:

> (x <- array(1:4, c(2,2)))
 [,1] [,2]
[1,]13
[2,]24

> as.array(x)
 [,1] [,2]
[1,]13
[2,]24

In version 2.8.0:
 
> (x <- array(1:4, c(2,2)))
 [,1] [,2]
[1,]13
[2,]24

> as.array(x)
Error: evaluation nested too deeply: infinite recursion /
options(expressions=)?


-
The power of accurate observation is commonly called cynicism 
by those who have not got it.  George Bernard Shaw
-- 
View this message in context: 
http://www.nabble.com/array-in-version-2.8.0-tp20295881p20295881.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] array in version 2.8.0

2008-11-02 Thread David Stoffer

Ok- it works now after flushing my objects.  The package Oarray doesn't seem
to work under version 2.8.0.  I emailed Jonathan about it and he said he
won't get to updating it until next year.  In the meantime, I was trying to
produce a seat-of-the-pants workaround and I probably redefined "as.array"
in my carelessness.  Sorry for waking everybody up.





Rolf Turner-3 wrote:
> 
> 
> On 3/11/2008, at 2:11 PM, David Stoffer wrote:
> 
>>
>> What happened?  TIA.
>>
>> In version 2.7.x:
>>
>>> (x <- array(1:4, c(2,2)))
>>  [,1] [,2]
>> [1,]13
>> [2,]24
>>
>>> as.array(x)
>>  [,1] [,2]
>> [1,]13
>> [2,]24
>>
>> In version 2.8.0:
>>
>>> (x <- array(1:4, c(2,2)))
>>  [,1] [,2]
>> [1,]13
>> [2,]24
>>
>>> as.array(x)
>> Error: evaluation nested too deeply: infinite recursion /
>> options(expressions=)?
> 
> Doesn't happen to me:
> 
>   > x <- array(1:4,c(2,2))
>   > x
>   [,1] [,2]
> [1,]13
> [2,]24
>   > class(x)
> [1] "matrix"
>   > as.array(x)
>   [,1] [,2]
> [1,]13
> [2,]24
> 
> My sessionInfo() is:
> 
>   > sessionInfo()
> R version 2.8.0 (2008-10-20)
> i386-apple-darwin8.11.1
> 
> locale:
> C
> 
> attached base packages:
> [1] datasets  utils stats graphics  grDevices methods   base
> 
> other attached packages:
> [1] misc_0.0-9 fortunes_1.3-5 MASS_7.2-44
> 
> Is there another array() function hanging around in your search path?
> 
>   cheers,
> 
>   Rolf Turner
> 
> 
> ##
> Attention:\ This e-mail message is privileged and confid...{{dropped:9}}
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 


-
The power of accurate observation is commonly called cynicism 
by those who have not got it.  George Bernard Shaw
-- 
View this message in context: 
http://www.nabble.com/array-in-version-2.8.0-tp20295881p20296616.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] "xreg" in ARIMA modelling.

2008-11-27 Thread David Stoffer

The help file states: "The exact likelihood is computed via a state-space
representation of the ARIMA process, and the innovations and their variance
found by a Kalman filter."   It is possible to include exogenous variables
(xreg) this way, but one can only assume this is done [only one person knows
for sure... the person who wrote the final version of arima(), and I hope he
chimes in to this].  If this is the case, then there is a likelihood
evaluation and AIC [or similar criteria, e.g., BIC and so on] would apply.


00alastair00 wrote:
> 
> Hello, 
> Does anyone know how the parameter estimates are calculated for xreg
> variables when called as part of an arima() command, or know of any
> literature that provides this info?  In particular, I was wondering if
> there is a quick way to compare different combinations of "xreg" variables
> in the arima() fit in the same way that you would in multiple regression
> (using AIC & R^2 etc.).  
> 
> Thanks very much!
> 
> 
> 


-
The power of accurate observation is commonly called cynicism 
by those who have not got it.  George Bernard Shaw
-- 
View this message in context: 
http://www.nabble.com/%22xreg%22-in-ARIMA-modelling.-tp20718058p20727625.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Spectral Analysis of Time Series in R

2008-12-05 Thread David Stoffer

You can do (1) and (2) [with some additional coding] using mvspec.R, which
you can download from http://www.stat.pitt.edu/stoffer/tsa2/chap7.htm ...
scroll down to the Spectral Envelope section and you'll find it there.  You
can look at the top part of the examples to get an idea of how to use
mvspec.R ... once you have the  spectral matrix estimate, you can code up
the extraction of the partial coherence and so on.





Alexander Schnebel wrote:
> 
> Dear R Community,
> 
> I am currently student at the Vienna University of Technology writing my 
> Diploma thesis on causality in time series and doing some analyses of 
> time series in R. I have the following questions:
> 
> (1) Is there a function in R to estimate the PARTIAL spectral coherence 
> of a multivariate time series? If yes, how does this work? Is there an 
> test in R if the partial spectral coherence between two variables is 
> zero? The functions I know (spectrum, etc.) only work to estimate the 
> spectral coherence.
> 
> (2) For some causality analysis I need an estimate of the inverse of the 
> spectral density matrix of a multivariate time series. Is there any 
> possibility in R to get this? Actually, I would be happy if I could at 
> least get a functional estimate of the spectral density matrix. I guess 
> this should work because R can plot the kernel density estimator of the 
> spectral density, so it should be possible to extract the underlying 
> function estimate.
> 
> (3) Is there any possibility to do Granger Causality in R? That means 
> fitting an VAR model and testing if some coefficients are zero.
> 
> Thank you very much in advance!
> 
> Best Regards,
> Alexander
> T
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 


-
The power of accurate observation is commonly called cynicism 
by those who have not got it.  George Bernard Shaw
-- 
View this message in context: 
http://www.nabble.com/Spectral-Analysis-of-Time-Series-in-R-tp20814256p20866634.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] arima, xreg, and the armax model

2009-05-03 Thread David Stoffer

Hi Marc- I have been [and am] extremely busy and haven't had much time to be
a playeR (lately I've become more of a moveR and shakeR ... some say more of
a boozeR and a loseR ... it's all prespective :).  I've updated the web page
with a little more info, but when I find the time I'll put up some ARMAX
examples on the Ch 6 page (state space models).  I think Paul Gilbert's DSE
package will do the job if you need it now.

The arima() help file says xreg does regression with autocorrelated errors,
and that is indeed what it does. The key line in arima.R is:  if(ncxreg > 0)
x <- x - xreg %*% par[narma + (1L:ncxreg)], so you see it replaces "x(t)"
with "x(t)-beta*xreg(t)" so to speak.  

... and, as usual, I stand on the GNU General Public License clause that
proclaims anything I say or do comes without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE (not that there's
anything wrong with being fit for a particular purpose).



Marc Vinyes-2 wrote:
> 
> Hello all,
>  
> I'm having fun again with the arima function. This time I read in:
> http://www.stat.pitt.edu/stoffer/tsa2/R_time_series_quick_fix.htm
>  
> < arima() does NOT fit an ARMAX model [insert slap head icon here]. This
> will
> be investigated as soon as time permits.>>
> (by R.H. Shumway & D.S. Stoffer)
>  
> This is quite surprising... Does anybody know anything about it?
>  
> Marc Vinyes (AleaSoft)
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 
> 


-
The power of accurate observation is commonly called cynicism 
by those who have not got it.  George Bernard Shaw
-- 
View this message in context: 
http://www.nabble.com/arima%2C-xreg%2C-and-the-armax-model-tp22724290p23361962.html
Sent from the R help mailing list archive at Nabble.com.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.