R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(1.79,194.9,1.95,195.5,2.26,196.0,2.04,196.2,2.16,196.2,2.75,196.2,2.79,196.2,2.88,197.0,3.36,197.7,2.97,198.0,3.10,198.2,2.49,198.5,2.2,198.6,2.25,199.5,2.09,200,2.79,201.3,3.14,202.2,2.93,202.9,2.65,203.5,2.67,203.5,2.26,204,2.35,204.1,2.13,204.3,2.18,204.5,2.9,204.8,2.63,205.1,2.67,205.7,1.81,206.5,1.33,206.9,0.88,207.1,1.28,207.8,1.26,208,1.26,208.5,1.29,208.6,1.1,209,1.37,209.1,1.21,209.7,1.74,209.8,1.76,209.9,1.48,210,1.04,210.8,1.62,211.4,1.49,211.7,1.79,212,1.8,212.2,1.58,212.4,1.86,212.9,1.74,213.4,1.59,213.7,1.26,214,1.13,214.3,1.92,214.8,2.61,215,2.26,215.9,2.41,216.4,2.26,216.9,2.03,217.2,2.86,217.5,2.55,217.9,2.27,218.1,2.26,218.6,2.57,218.9,3.07,219.3,2.76,220.4,2.51,220.9,2.87,221,3.14,221.8,3.11,222,3.16,222.2,2.47,222.5,2.57,222.9,2.89,223.1),dim=c(2,72),dimnames=list(c('Yt','Xt'),1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('Yt','Xt'),1:72))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Yt Xt M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 1.79 194.9 1 0 0 0 0 0 0 0 0 0 0
2 1.95 195.5 0 1 0 0 0 0 0 0 0 0 0
3 2.26 196.0 0 0 1 0 0 0 0 0 0 0 0
4 2.04 196.2 0 0 0 1 0 0 0 0 0 0 0
5 2.16 196.2 0 0 0 0 1 0 0 0 0 0 0
6 2.75 196.2 0 0 0 0 0 1 0 0 0 0 0
7 2.79 196.2 0 0 0 0 0 0 1 0 0 0 0
8 2.88 197.0 0 0 0 0 0 0 0 1 0 0 0
9 3.36 197.7 0 0 0 0 0 0 0 0 1 0 0
10 2.97 198.0 0 0 0 0 0 0 0 0 0 1 0
11 3.10 198.2 0 0 0 0 0 0 0 0 0 0 1
12 2.49 198.5 0 0 0 0 0 0 0 0 0 0 0
13 2.20 198.6 1 0 0 0 0 0 0 0 0 0 0
14 2.25 199.5 0 1 0 0 0 0 0 0 0 0 0
15 2.09 200.0 0 0 1 0 0 0 0 0 0 0 0
16 2.79 201.3 0 0 0 1 0 0 0 0 0 0 0
17 3.14 202.2 0 0 0 0 1 0 0 0 0 0 0
18 2.93 202.9 0 0 0 0 0 1 0 0 0 0 0
19 2.65 203.5 0 0 0 0 0 0 1 0 0 0 0
20 2.67 203.5 0 0 0 0 0 0 0 1 0 0 0
21 2.26 204.0 0 0 0 0 0 0 0 0 1 0 0
22 2.35 204.1 0 0 0 0 0 0 0 0 0 1 0
23 2.13 204.3 0 0 0 0 0 0 0 0 0 0 1
24 2.18 204.5 0 0 0 0 0 0 0 0 0 0 0
25 2.90 204.8 1 0 0 0 0 0 0 0 0 0 0
26 2.63 205.1 0 1 0 0 0 0 0 0 0 0 0
27 2.67 205.7 0 0 1 0 0 0 0 0 0 0 0
28 1.81 206.5 0 0 0 1 0 0 0 0 0 0 0
29 1.33 206.9 0 0 0 0 1 0 0 0 0 0 0
30 0.88 207.1 0 0 0 0 0 1 0 0 0 0 0
31 1.28 207.8 0 0 0 0 0 0 1 0 0 0 0
32 1.26 208.0 0 0 0 0 0 0 0 1 0 0 0
33 1.26 208.5 0 0 0 0 0 0 0 0 1 0 0
34 1.29 208.6 0 0 0 0 0 0 0 0 0 1 0
35 1.10 209.0 0 0 0 0 0 0 0 0 0 0 1
36 1.37 209.1 0 0 0 0 0 0 0 0 0 0 0
37 1.21 209.7 1 0 0 0 0 0 0 0 0 0 0
38 1.74 209.8 0 1 0 0 0 0 0 0 0 0 0
39 1.76 209.9 0 0 1 0 0 0 0 0 0 0 0
40 1.48 210.0 0 0 0 1 0 0 0 0 0 0 0
41 1.04 210.8 0 0 0 0 1 0 0 0 0 0 0
42 1.62 211.4 0 0 0 0 0 1 0 0 0 0 0
43 1.49 211.7 0 0 0 0 0 0 1 0 0 0 0
44 1.79 212.0 0 0 0 0 0 0 0 1 0 0 0
45 1.80 212.2 0 0 0 0 0 0 0 0 1 0 0
46 1.58 212.4 0 0 0 0 0 0 0 0 0 1 0
47 1.86 212.9 0 0 0 0 0 0 0 0 0 0 1
48 1.74 213.4 0 0 0 0 0 0 0 0 0 0 0
49 1.59 213.7 1 0 0 0 0 0 0 0 0 0 0
50 1.26 214.0 0 1 0 0 0 0 0 0 0 0 0
51 1.13 214.3 0 0 1 0 0 0 0 0 0 0 0
52 1.92 214.8 0 0 0 1 0 0 0 0 0 0 0
53 2.61 215.0 0 0 0 0 1 0 0 0 0 0 0
54 2.26 215.9 0 0 0 0 0 1 0 0 0 0 0
55 2.41 216.4 0 0 0 0 0 0 1 0 0 0 0
56 2.26 216.9 0 0 0 0 0 0 0 1 0 0 0
57 2.03 217.2 0 0 0 0 0 0 0 0 1 0 0
58 2.86 217.5 0 0 0 0 0 0 0 0 0 1 0
59 2.55 217.9 0 0 0 0 0 0 0 0 0 0 1
60 2.27 218.1 0 0 0 0 0 0 0 0 0 0 0
61 2.26 218.6 1 0 0 0 0 0 0 0 0 0 0
62 2.57 218.9 0 1 0 0 0 0 0 0 0 0 0
63 3.07 219.3 0 0 1 0 0 0 0 0 0 0 0
64 2.76 220.4 0 0 0 1 0 0 0 0 0 0 0
65 2.51 220.9 0 0 0 0 1 0 0 0 0 0 0
66 2.87 221.0 0 0 0 0 0 1 0 0 0 0 0
67 3.14 221.8 0 0 0 0 0 0 1 0 0 0 0
68 3.11 222.0 0 0 0 0 0 0 0 1 0 0 0
69 3.16 222.2 0 0 0 0 0 0 0 0 1 0 0
70 2.47 222.5 0 0 0 0 0 0 0 0 0 1 0
71 2.57 222.9 0 0 0 0 0 0 0 0 0 0 1
72 2.89 223.1 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Xt M1 M2 M3 M4
2.2929372 -0.0006455 -0.1678401 -0.0925711 0.0043537 -0.0252160
M5 M6 M7 M8 M9 M10
-0.0265814 0.0603542 0.1356662 0.1708813 0.1544729 0.0962794
M11
0.0615053
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.33961 -0.43900 0.06764 0.51403 1.04020
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.2929372 2.1228813 1.080 0.284
Xt -0.0006455 0.0099664 -0.065 0.949
M1 -0.1678401 0.4010866 -0.418 0.677
M2 -0.0925711 0.4006539 -0.231 0.818
M3 0.0043537 0.4002785 0.011 0.991
M4 -0.0252160 0.3997405 -0.063 0.950
M5 -0.0265814 0.3994292 -0.067 0.947
M6 0.0603542 0.3991969 0.151 0.880
M7 0.1356662 0.3989813 0.340 0.735
M8 0.1708813 0.3988665 0.428 0.670
M9 0.1544729 0.3987652 0.387 0.700
M10 0.0962794 0.3987270 0.241 0.810
M11 0.0615053 0.3986899 0.154 0.878
Residual standard error: 0.6905 on 59 degrees of freedom
Multiple R-squared: 0.02364, Adjusted R-squared: -0.1749
F-statistic: 0.119 on 12 and 59 DF, p-value: 0.9999
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.05648573 0.112971459 0.943514270
[2,] 0.03799457 0.075989143 0.962005429
[3,] 0.03246627 0.064932547 0.967533727
[4,] 0.04337938 0.086758761 0.956620619
[5,] 0.04227744 0.084554876 0.957722562
[6,] 0.17233043 0.344660856 0.827669572
[7,] 0.17845825 0.356916506 0.821541747
[8,] 0.23336021 0.466720419 0.766639791
[9,] 0.20009841 0.400196820 0.799901590
[10,] 0.51729250 0.965415002 0.482707501
[11,] 0.71987252 0.560254966 0.280127483
[12,] 0.91945078 0.161098442 0.080549221
[13,] 0.95658749 0.086825025 0.043412512
[14,] 0.98752622 0.024947554 0.012473777
[15,] 0.99857049 0.002859018 0.001429509
[16,] 0.99884711 0.002305779 0.001152890
[17,] 0.99894980 0.002100406 0.001050203
[18,] 0.99889839 0.002203216 0.001101608
[19,] 0.99846104 0.003077918 0.001538959
[20,] 0.99801142 0.003977166 0.001988583
[21,] 0.99646293 0.007074148 0.003537074
[22,] 0.99351235 0.012975291 0.006487645
[23,] 0.99343383 0.013132338 0.006566169
[24,] 0.99224210 0.015515794 0.007757897
[25,] 0.98636761 0.027264773 0.013632387
[26,] 0.98559004 0.028819917 0.014409958
[27,] 0.97516735 0.049665304 0.024832652
[28,] 0.96406394 0.071872125 0.035936063
[29,] 0.94214541 0.115709186 0.057854593
[30,] 0.91121852 0.177562963 0.088781481
[31,] 0.86747536 0.265049290 0.132524645
[32,] 0.82208478 0.355830440 0.177915220
[33,] 0.75679737 0.486405251 0.243202626
[34,] 0.67632385 0.647352291 0.323676145
[35,] 0.65098021 0.698039584 0.349019792
[36,] 0.88190169 0.236196622 0.118098311
[37,] 0.84267615 0.314647698 0.157323849
[38,] 0.87348001 0.253039973 0.126519986
[39,] 0.79887957 0.402240857 0.201120429
[40,] 0.70325963 0.593480743 0.296740371
[41,] 0.59953012 0.800939760 0.400469880
> postscript(file="/var/www/html/rcomp/tmp/1jgk61291467829.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/2t71r1291467829.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/3t71r1291467829.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/4t71r1291467829.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/5mg0b1291467829.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 72
Frequency = 1
1 2 3 4 5 6
-0.20929403 -0.12417569 0.08922219 -0.10107903 0.02028641 0.52335080
7 8 9 10 11 12
0.48803882 0.54334004 1.04020035 0.70858747 0.87349065 0.32518959
13 14 15 16 17 18
0.20309423 0.17840621 -0.07819591 0.65221289 1.00415926 0.70767548
19 20 21 22 23 24
0.35275079 0.33753563 -0.05573316 0.09252487 -0.09257195 0.01906244
25 26 27 28 29 30
0.90709617 0.56202087 0.50548330 -0.32443064 -0.80280701 -1.33961353
31 32 33 34 35 36
-1.01447367 -1.06955974 -1.05282852 -0.96457049 -1.11953822 -0.78796837
37 38 39 40 41 42
-0.77974100 -0.32494540 -0.40180571 -0.65217148 -1.09028965 -0.59683798
43 44 45 46 47 48
-0.80195632 -0.53697784 -0.51044026 -0.67211769 -0.35702087 -0.41519283
49 50 51 52 53 54
-0.39715910 -0.80223441 -1.02896562 -0.20907320 0.48242134 0.04606665
55 56 57 58 59 60
0.12107741 -0.06381501 -0.27721289 0.61117423 0.33620651 0.11784090
61 62 63 64 65 66
0.27600373 0.51092842 0.91426175 0.63454146 0.38622964 0.65935858
67 68 69 70 71 72
0.85456298 0.78947691 0.85601448 0.22440161 0.35943388 0.74106827
> postscript(file="/var/www/html/rcomp/tmp/6mg0b1291467829.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.20929403 NA
1 -0.12417569 -0.20929403
2 0.08922219 -0.12417569
3 -0.10107903 0.08922219
4 0.02028641 -0.10107903
5 0.52335080 0.02028641
6 0.48803882 0.52335080
7 0.54334004 0.48803882
8 1.04020035 0.54334004
9 0.70858747 1.04020035
10 0.87349065 0.70858747
11 0.32518959 0.87349065
12 0.20309423 0.32518959
13 0.17840621 0.20309423
14 -0.07819591 0.17840621
15 0.65221289 -0.07819591
16 1.00415926 0.65221289
17 0.70767548 1.00415926
18 0.35275079 0.70767548
19 0.33753563 0.35275079
20 -0.05573316 0.33753563
21 0.09252487 -0.05573316
22 -0.09257195 0.09252487
23 0.01906244 -0.09257195
24 0.90709617 0.01906244
25 0.56202087 0.90709617
26 0.50548330 0.56202087
27 -0.32443064 0.50548330
28 -0.80280701 -0.32443064
29 -1.33961353 -0.80280701
30 -1.01447367 -1.33961353
31 -1.06955974 -1.01447367
32 -1.05282852 -1.06955974
33 -0.96457049 -1.05282852
34 -1.11953822 -0.96457049
35 -0.78796837 -1.11953822
36 -0.77974100 -0.78796837
37 -0.32494540 -0.77974100
38 -0.40180571 -0.32494540
39 -0.65217148 -0.40180571
40 -1.09028965 -0.65217148
41 -0.59683798 -1.09028965
42 -0.80195632 -0.59683798
43 -0.53697784 -0.80195632
44 -0.51044026 -0.53697784
45 -0.67211769 -0.51044026
46 -0.35702087 -0.67211769
47 -0.41519283 -0.35702087
48 -0.39715910 -0.41519283
49 -0.80223441 -0.39715910
50 -1.02896562 -0.80223441
51 -0.20907320 -1.02896562
52 0.48242134 -0.20907320
53 0.04606665 0.48242134
54 0.12107741 0.04606665
55 -0.06381501 0.12107741
56 -0.27721289 -0.06381501
57 0.61117423 -0.27721289
58 0.33620651 0.61117423
59 0.11784090 0.33620651
60 0.27600373 0.11784090
61 0.51092842 0.27600373
62 0.91426175 0.51092842
63 0.63454146 0.91426175
64 0.38622964 0.63454146
65 0.65935858 0.38622964
66 0.85456298 0.65935858
67 0.78947691 0.85456298
68 0.85601448 0.78947691
69 0.22440161 0.85601448
70 0.35943388 0.22440161
71 0.74106827 0.35943388
72 NA 0.74106827
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.12417569 -0.20929403
[2,] 0.08922219 -0.12417569
[3,] -0.10107903 0.08922219
[4,] 0.02028641 -0.10107903
[5,] 0.52335080 0.02028641
[6,] 0.48803882 0.52335080
[7,] 0.54334004 0.48803882
[8,] 1.04020035 0.54334004
[9,] 0.70858747 1.04020035
[10,] 0.87349065 0.70858747
[11,] 0.32518959 0.87349065
[12,] 0.20309423 0.32518959
[13,] 0.17840621 0.20309423
[14,] -0.07819591 0.17840621
[15,] 0.65221289 -0.07819591
[16,] 1.00415926 0.65221289
[17,] 0.70767548 1.00415926
[18,] 0.35275079 0.70767548
[19,] 0.33753563 0.35275079
[20,] -0.05573316 0.33753563
[21,] 0.09252487 -0.05573316
[22,] -0.09257195 0.09252487
[23,] 0.01906244 -0.09257195
[24,] 0.90709617 0.01906244
[25,] 0.56202087 0.90709617
[26,] 0.50548330 0.56202087
[27,] -0.32443064 0.50548330
[28,] -0.80280701 -0.32443064
[29,] -1.33961353 -0.80280701
[30,] -1.01447367 -1.33961353
[31,] -1.06955974 -1.01447367
[32,] -1.05282852 -1.06955974
[33,] -0.96457049 -1.05282852
[34,] -1.11953822 -0.96457049
[35,] -0.78796837 -1.11953822
[36,] -0.77974100 -0.78796837
[37,] -0.32494540 -0.77974100
[38,] -0.40180571 -0.32494540
[39,] -0.65217148 -0.40180571
[40,] -1.09028965 -0.65217148
[41,] -0.59683798 -1.09028965
[42,] -0.80195632 -0.59683798
[43,] -0.53697784 -0.80195632
[44,] -0.51044026 -0.53697784
[45,] -0.67211769 -0.51044026
[46,] -0.35702087 -0.67211769
[47,] -0.41519283 -0.35702087
[48,] -0.39715910 -0.41519283
[49,] -0.80223441 -0.39715910
[50,] -1.02896562 -0.80223441
[51,] -0.20907320 -1.02896562
[52,] 0.48242134 -0.20907320
[53,] 0.04606665 0.48242134
[54,] 0.12107741 0.04606665
[55,] -0.06381501 0.12107741
[56,] -0.27721289 -0.06381501
[57,] 0.61117423 -0.27721289
[58,] 0.33620651 0.61117423
[59,] 0.11784090 0.33620651
[60,] 0.27600373 0.11784090
[61,] 0.51092842 0.27600373
[62,] 0.91426175 0.51092842
[63,] 0.63454146 0.91426175
[64,] 0.38622964 0.63454146
[65,] 0.65935858 0.38622964
[66,] 0.85456298 0.65935858
[67,] 0.78947691 0.85456298
[68,] 0.85601448 0.78947691
[69,] 0.22440161 0.85601448
[70,] 0.35943388 0.22440161
[71,] 0.74106827 0.35943388
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.12417569 -0.20929403
2 0.08922219 -0.12417569
3 -0.10107903 0.08922219
4 0.02028641 -0.10107903
5 0.52335080 0.02028641
6 0.48803882 0.52335080
7 0.54334004 0.48803882
8 1.04020035 0.54334004
9 0.70858747 1.04020035
10 0.87349065 0.70858747
11 0.32518959 0.87349065
12 0.20309423 0.32518959
13 0.17840621 0.20309423
14 -0.07819591 0.17840621
15 0.65221289 -0.07819591
16 1.00415926 0.65221289
17 0.70767548 1.00415926
18 0.35275079 0.70767548
19 0.33753563 0.35275079
20 -0.05573316 0.33753563
21 0.09252487 -0.05573316
22 -0.09257195 0.09252487
23 0.01906244 -0.09257195
24 0.90709617 0.01906244
25 0.56202087 0.90709617
26 0.50548330 0.56202087
27 -0.32443064 0.50548330
28 -0.80280701 -0.32443064
29 -1.33961353 -0.80280701
30 -1.01447367 -1.33961353
31 -1.06955974 -1.01447367
32 -1.05282852 -1.06955974
33 -0.96457049 -1.05282852
34 -1.11953822 -0.96457049
35 -0.78796837 -1.11953822
36 -0.77974100 -0.78796837
37 -0.32494540 -0.77974100
38 -0.40180571 -0.32494540
39 -0.65217148 -0.40180571
40 -1.09028965 -0.65217148
41 -0.59683798 -1.09028965
42 -0.80195632 -0.59683798
43 -0.53697784 -0.80195632
44 -0.51044026 -0.53697784
45 -0.67211769 -0.51044026
46 -0.35702087 -0.67211769
47 -0.41519283 -0.35702087
48 -0.39715910 -0.41519283
49 -0.80223441 -0.39715910
50 -1.02896562 -0.80223441
51 -0.20907320 -1.02896562
52 0.48242134 -0.20907320
53 0.04606665 0.48242134
54 0.12107741 0.04606665
55 -0.06381501 0.12107741
56 -0.27721289 -0.06381501
57 0.61117423 -0.27721289
58 0.33620651 0.61117423
59 0.11784090 0.33620651
60 0.27600373 0.11784090
61 0.51092842 0.27600373
62 0.91426175 0.51092842
63 0.63454146 0.91426175
64 0.38622964 0.63454146
65 0.65935858 0.38622964
66 0.85456298 0.65935858
67 0.78947691 0.85456298
68 0.85601448 0.78947691
69 0.22440161 0.85601448
70 0.35943388 0.22440161
71 0.74106827 0.35943388
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/7fq0x1291467829.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/8fq0x1291467829.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/9fq0x1291467829.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/rcomp/tmp/107zh01291467829.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/11b0fo1291467829.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/12w0wt1291467829.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/13sau21291467829.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/14ess81291467829.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/15htrw1291467829.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/162t7k1291467829.tab")
+ }
>
> try(system("convert tmp/1jgk61291467829.ps tmp/1jgk61291467829.png",intern=TRUE))
character(0)
> try(system("convert tmp/2t71r1291467829.ps tmp/2t71r1291467829.png",intern=TRUE))
character(0)
> try(system("convert tmp/3t71r1291467829.ps tmp/3t71r1291467829.png",intern=TRUE))
character(0)
> try(system("convert tmp/4t71r1291467829.ps tmp/4t71r1291467829.png",intern=TRUE))
character(0)
> try(system("convert tmp/5mg0b1291467829.ps tmp/5mg0b1291467829.png",intern=TRUE))
character(0)
> try(system("convert tmp/6mg0b1291467829.ps tmp/6mg0b1291467829.png",intern=TRUE))
character(0)
> try(system("convert tmp/7fq0x1291467829.ps tmp/7fq0x1291467829.png",intern=TRUE))
character(0)
> try(system("convert tmp/8fq0x1291467829.ps tmp/8fq0x1291467829.png",intern=TRUE))
character(0)
> try(system("convert tmp/9fq0x1291467829.ps tmp/9fq0x1291467829.png",intern=TRUE))
character(0)
> try(system("convert tmp/107zh01291467829.ps tmp/107zh01291467829.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.606 1.645 8.264