## Thursday, December 27, 2018

### Estimating Linear Trend -R

Estimating a Linear Trend
x = arima.sim(list(order=c(2,0,0), ar=c(1,-.9)), n=2^8)
> (u = polyroot(c(1,-1,.9)))
[1] 0.5555556+0.8958064i 0.5555556-0.8958064i
> Arg(u[1])/(2*pi)
[1] 0.1616497
> par(mfcol=c(4,1))
> plot.ts(x)
> spec.pgram(x,spans=c(3,3),log="no")
> spec.ar(x, log="no")
> spec.arma(ar=c(1,-.9), log="no")
> install.packages("spec.arma")
Warning in install.packages
> spec.pgram(x,
+            taper=0, fast=FALSE, detrend=FALSE, log="no")
> spec.arma(ar=c(1,-.9), log="no")
> per = abs(fft(x))^2/length(x)
> set.seed(666)
> x = 50 + arima.sim(list(order=c(1,0,1), ar=.9, ma=-.5), n=200)
> acf(x); pacf(x)
> (x.fit = arima(x, order = c(1, 0, 1)))

Call:
arima(x = x, order = c(1, 0, 1))

Coefficients:
ar1     ma1  intercept
0.8340  -0.432    49.8962
s.e.  0.0645   0.111     0.2452

sigma^2 estimated as 1.07:  log likelihood = -290.79,  aic = 589.58
> tsdiag(x.fit, gof.lag=20)
> (x.ima = HoltWinters(x, beta=FALSE, gamma=FALSE))
Holt-Winters exponential smoothing without trend and without seasonal component.

Call:
HoltWinters(x = x, beta = FALSE, gamma = FALSE)

Smoothing parameters:
alpha: 0.4240096
beta : FALSE
gamma: FALSE

Coefficients:
[,1]
a 49.45997
> plot(x.ima)
> set.seed(111)
> phi.yw = rep(NA, 1000)
> for (i in 1:1000){
+      e = rexp(150, rate=.5); u = runif(150,-1,1); de = e*sign(u)
+      x = 50 + arima.sim(n=100,list(ar=.95), innov=de, n.start=50)
+      phi.yw[i] = ar.yw(x, order=1)\$ar }
> hist(phi.yw, prob=TRUE, main="")
> lines(density(phi.yw, bw=.015))
> ACF = ARMAacf(ar=c(1.5,-.75), ma=0, 24)[-1]
> PACF = ARMAacf(ar=c(1.5,-.75), ma=0, 24, pacf=TRUE)
> par(mfrow=c(1,2))
> plot(ACF, type="h", xlab="lag", ylim=c(-.8,1)); abline(h=0)
> plot(PACF, type="h", xlab="lag", ylim=c(-.8,1)); abline(h=0)
> ARMAtoMA(ar=.9, ma=.5, 50)
[1] 1.400000000 1.260000000 1.134000000
[4] 1.020600000 0.918540000 0.826686000
[7] 0.744017400 0.669615660 0.602654094
[10] 0.542388685 0.488149816 0.439334835
[13] 0.395401351 0.355861216 0.320275094
[16] 0.288247585 0.259422826 0.233480544
[19] 0.210132489 0.189119240 0.170207316
[22] 0.153186585 0.137867926 0.124081134
[25] 0.111673020 0.100505718 0.090455146
[28] 0.081409632 0.073268669 0.065941802
[31] 0.059347622 0.053412859 0.048071573
[34] 0.043264416 0.038937975 0.035044177
[37] 0.031539759 0.028385783 0.025547205
[40] 0.022992485 0.020693236 0.018623913
[43] 0.016761521 0.015085369 0.013576832
[46] 0.012219149 0.010997234 0.009897511
[49] 0.008907760 0.008016984
> plot(ARMAtoMA(ar=.9, ma=.5, 50))
> z = c(1,-1.5,.75)
> (a = polyroot(z)[1])
[1] 1+0.57735i
> arg = Arg(a)/(2*pi)
> 1/arg
[1] 12
> set.seed(90210)
> ar2 = arima.sim(list(order=c(2,0,0), ar=c(1.5,-.75)), n = 144)
> plot(1:144/12, ar2, type="l", xlab="Time (one unit = 12 points)")
> abline(v=0:12, lty="dotted", lwd=2)
> ACF = ARMAacf(ar=c(1.5,-.75), ma=0, 50)
> plot(ACF, type="h", xlab="lag")
> abline(h=0)
> par(mfrow=c(2,1)
+ plot(arima.sim(list(order=c(1,0,0), ar=.9), n=100), ylab="x",
Error: unexpected symbol in:
"par(mfrow=c(2,1)
plot"
>      main=(expression(AR(1)~~~phi==+.9))
> plot(arima.sim(list(order=c(1,0,0), ar=.9), n=100), ylab="x",main=(expression(AR(1)~~~phi==+.9))
+
+ r
Error: unexpected symbol in:
"
r"
> plot(arima.sim(list(order=c(1,0,0), ar=-.9), n=100), ylab="x",
+      main=(expression(AR(1)~~~phi==-.9))
+
+ b
Error: unexpected symbol in:
"
b"
> par(mfrow = c(2,1))
> plot(arima.sim(list(order=c(0,0,1), ma=.5), n=100), ylab="x",
+      main=(expression(MA(1)~~~theta==+.5)))
> plot(arima.sim(list(order=c(0,0,1), ma=-.5), n=100), ylab="x",
+      main=(expression(MA(1)~~~theta==-.5)))
> par(mfrow=c(2,1))
> set.seed(1000)
> x = 2*cos(2*pi*1:500/50 + .6*pi) + rnorm(500,0,5)
> z1 = cos(2*pi*1:500/50); z2 = sin(2*pi*1:500/50)
> summary(fit <- lm(x~0+z1+z2))

Call:
lm(formula = x ~ 0 + z1 + z2)

Residuals:
Min       1Q   Median       3Q
-12.6589  -3.0872   0.1508   3.1302
Max
13.5821

Coefficients:
Estimate Std. Error t value Pr(>|t|)
z1  -0.7083     0.2958  -2.394    0.017
z2  -2.5473     0.2958  -8.611   <2e-16

z1 *
z2 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’
0.1 ‘ ’ 1

Residual standard error: 4.678 on 498 degrees of freedom
Multiple R-squared:  0.1382, Adjusted R-squared:  0.1348
F-statistic: 39.94 on 2 and 498 DF,  p-value: < 2.2e-16

> plot.ts(x, lty="dashed")
> lines(fitted(fit), lwd=2)
> I = abs(fft(x))^2/500
> P = (4/500)*I[1:250]
> f = 0:249/500
> plot(f, P, type="l", xlab="Frequency", ylab="Scaled Periodogram")