R version 2.7.0 (2008-04-22)
Copyright (C) 2008 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(2648.9,0,2669.6,0,3042.3,0,2604.2,0,2732.1,0,2621.7,0,2483.7,0,2479.3,0,2684.6,0,2834.7,0,2566.1,0,2251.2,0,2350,1,2299.8,1,2542.8,1,2530.2,1,2508.1,1,2616.8,1,2534.1,1,2181.8,1,2578.9,1,2841.9,1,2529.9,1,2103.2,1,2326.2,1,2452.6,1,2782.1,1,2727.3,1,2648.2,1,2760.7,1,2613,1,2225.4,1,2713.9,1,2923.3,1,2707,1,2473.9,1,2521,1,2531.8,1,3068.8,1,2826.9,1,2674.2,1,2966.6,1,2798.8,1,2629.6,1,3124.6,1,3115.7,1,3083,1,2863.9,1,2728.7,1,2789.4,1,3225.7,1,3148.2,1,2836.5,1,3153.5,1,2656.9,1,2834.7,1,3172.5,1,2998.8,1,3103.1,1,2735.6,1,2818.1,1,2874.4,1,3438.5,1,2949.1,1,3306.8,1,3530,1,3003.8,1,3206.4,1,3514.6,1,3522.6,1,3525.5,1,2996.2,1,3231.1,1,3030,1,3541.7,1,3113.2,1,3390.8,1,3424.2,1,3079.8,1,3123.4,1,3317.1,1,3579.9,1,3317.9,1,2668.1,1),dim=c(2,84),dimnames=list(c('Y','X'),1:84))
> y <- array(NA,dim=c(2,84),dimnames=list(c('Y','X'),1:84))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = '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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 2648.9 0 1 0 0 0 0 0 0 0 0 0 0 1
2 2669.6 0 0 1 0 0 0 0 0 0 0 0 0 2
3 3042.3 0 0 0 1 0 0 0 0 0 0 0 0 3
4 2604.2 0 0 0 0 1 0 0 0 0 0 0 0 4
5 2732.1 0 0 0 0 0 1 0 0 0 0 0 0 5
6 2621.7 0 0 0 0 0 0 1 0 0 0 0 0 6
7 2483.7 0 0 0 0 0 0 0 1 0 0 0 0 7
8 2479.3 0 0 0 0 0 0 0 0 1 0 0 0 8
9 2684.6 0 0 0 0 0 0 0 0 0 1 0 0 9
10 2834.7 0 0 0 0 0 0 0 0 0 0 1 0 10
11 2566.1 0 0 0 0 0 0 0 0 0 0 0 1 11
12 2251.2 0 0 0 0 0 0 0 0 0 0 0 0 12
13 2350.0 1 1 0 0 0 0 0 0 0 0 0 0 13
14 2299.8 1 0 1 0 0 0 0 0 0 0 0 0 14
15 2542.8 1 0 0 1 0 0 0 0 0 0 0 0 15
16 2530.2 1 0 0 0 1 0 0 0 0 0 0 0 16
17 2508.1 1 0 0 0 0 1 0 0 0 0 0 0 17
18 2616.8 1 0 0 0 0 0 1 0 0 0 0 0 18
19 2534.1 1 0 0 0 0 0 0 1 0 0 0 0 19
20 2181.8 1 0 0 0 0 0 0 0 1 0 0 0 20
21 2578.9 1 0 0 0 0 0 0 0 0 1 0 0 21
22 2841.9 1 0 0 0 0 0 0 0 0 0 1 0 22
23 2529.9 1 0 0 0 0 0 0 0 0 0 0 1 23
24 2103.2 1 0 0 0 0 0 0 0 0 0 0 0 24
25 2326.2 1 1 0 0 0 0 0 0 0 0 0 0 25
26 2452.6 1 0 1 0 0 0 0 0 0 0 0 0 26
27 2782.1 1 0 0 1 0 0 0 0 0 0 0 0 27
28 2727.3 1 0 0 0 1 0 0 0 0 0 0 0 28
29 2648.2 1 0 0 0 0 1 0 0 0 0 0 0 29
30 2760.7 1 0 0 0 0 0 1 0 0 0 0 0 30
31 2613.0 1 0 0 0 0 0 0 1 0 0 0 0 31
32 2225.4 1 0 0 0 0 0 0 0 1 0 0 0 32
33 2713.9 1 0 0 0 0 0 0 0 0 1 0 0 33
34 2923.3 1 0 0 0 0 0 0 0 0 0 1 0 34
35 2707.0 1 0 0 0 0 0 0 0 0 0 0 1 35
36 2473.9 1 0 0 0 0 0 0 0 0 0 0 0 36
37 2521.0 1 1 0 0 0 0 0 0 0 0 0 0 37
38 2531.8 1 0 1 0 0 0 0 0 0 0 0 0 38
39 3068.8 1 0 0 1 0 0 0 0 0 0 0 0 39
40 2826.9 1 0 0 0 1 0 0 0 0 0 0 0 40
41 2674.2 1 0 0 0 0 1 0 0 0 0 0 0 41
42 2966.6 1 0 0 0 0 0 1 0 0 0 0 0 42
43 2798.8 1 0 0 0 0 0 0 1 0 0 0 0 43
44 2629.6 1 0 0 0 0 0 0 0 1 0 0 0 44
45 3124.6 1 0 0 0 0 0 0 0 0 1 0 0 45
46 3115.7 1 0 0 0 0 0 0 0 0 0 1 0 46
47 3083.0 1 0 0 0 0 0 0 0 0 0 0 1 47
48 2863.9 1 0 0 0 0 0 0 0 0 0 0 0 48
49 2728.7 1 1 0 0 0 0 0 0 0 0 0 0 49
50 2789.4 1 0 1 0 0 0 0 0 0 0 0 0 50
51 3225.7 1 0 0 1 0 0 0 0 0 0 0 0 51
52 3148.2 1 0 0 0 1 0 0 0 0 0 0 0 52
53 2836.5 1 0 0 0 0 1 0 0 0 0 0 0 53
54 3153.5 1 0 0 0 0 0 1 0 0 0 0 0 54
55 2656.9 1 0 0 0 0 0 0 1 0 0 0 0 55
56 2834.7 1 0 0 0 0 0 0 0 1 0 0 0 56
57 3172.5 1 0 0 0 0 0 0 0 0 1 0 0 57
58 2998.8 1 0 0 0 0 0 0 0 0 0 1 0 58
59 3103.1 1 0 0 0 0 0 0 0 0 0 0 1 59
60 2735.6 1 0 0 0 0 0 0 0 0 0 0 0 60
61 2818.1 1 1 0 0 0 0 0 0 0 0 0 0 61
62 2874.4 1 0 1 0 0 0 0 0 0 0 0 0 62
63 3438.5 1 0 0 1 0 0 0 0 0 0 0 0 63
64 2949.1 1 0 0 0 1 0 0 0 0 0 0 0 64
65 3306.8 1 0 0 0 0 1 0 0 0 0 0 0 65
66 3530.0 1 0 0 0 0 0 1 0 0 0 0 0 66
67 3003.8 1 0 0 0 0 0 0 1 0 0 0 0 67
68 3206.4 1 0 0 0 0 0 0 0 1 0 0 0 68
69 3514.6 1 0 0 0 0 0 0 0 0 1 0 0 69
70 3522.6 1 0 0 0 0 0 0 0 0 0 1 0 70
71 3525.5 1 0 0 0 0 0 0 0 0 0 0 1 71
72 2996.2 1 0 0 0 0 0 0 0 0 0 0 0 72
73 3231.1 1 1 0 0 0 0 0 0 0 0 0 0 73
74 3030.0 1 0 1 0 0 0 0 0 0 0 0 0 74
75 3541.7 1 0 0 1 0 0 0 0 0 0 0 0 75
76 3113.2 1 0 0 0 1 0 0 0 0 0 0 0 76
77 3390.8 1 0 0 0 0 1 0 0 0 0 0 0 77
78 3424.2 1 0 0 0 0 0 1 0 0 0 0 0 78
79 3079.8 1 0 0 0 0 0 0 1 0 0 0 0 79
80 3123.4 1 0 0 0 0 0 0 0 1 0 0 0 80
81 3317.1 1 0 0 0 0 0 0 0 0 1 0 0 81
82 3579.9 1 0 0 0 0 0 0 0 0 0 1 0 82
83 3317.9 1 0 0 0 0 0 0 0 0 0 0 1 83
84 2668.1 1 0 0 0 0 0 0 0 0 0 0 0 84
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
2201.40 -321.56 226.96 216.61 630.64 367.94
M5 M6 M7 M8 M9 M10
382.45 508.26 222.63 138.97 471.76 559.56
M11 t
405.21 13.73
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-364.652 -94.966 6.635 80.981 325.251
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2201.4038 66.4662 33.121 < 2e-16 ***
X -321.5589 56.0112 -5.741 2.23e-07 ***
M1 226.9616 76.4395 2.969 0.00409 **
M2 216.6079 76.3478 2.837 0.00595 **
M3 630.6400 76.2648 8.269 5.84e-12 ***
M4 367.9435 76.1904 4.829 7.79e-06 ***
M5 382.4470 76.1247 5.024 3.72e-06 ***
M6 508.2648 76.0678 6.682 4.73e-09 ***
M7 222.6254 76.0195 2.929 0.00459 **
M8 138.9717 75.9800 1.829 0.07165 .
M9 471.7610 75.9493 6.212 3.30e-08 ***
M10 559.5644 75.9274 7.370 2.63e-10 ***
M11 405.2108 75.9142 5.338 1.10e-06 ***
t 13.7251 0.8167 16.806 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 142 on 70 degrees of freedom
Multiple R-squared: 0.8652, Adjusted R-squared: 0.8402
F-statistic: 34.56 on 13 and 70 DF, p-value: < 2.2e-16
> 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.50681129 0.98637741 0.4931887
[2,] 0.58794390 0.82411219 0.4120561
[3,] 0.66565973 0.66868054 0.3343403
[4,] 0.56521410 0.86957180 0.4347859
[5,] 0.47270593 0.94541185 0.5272941
[6,] 0.45898855 0.91797710 0.5410114
[7,] 0.40166400 0.80332801 0.5983360
[8,] 0.31652243 0.63304486 0.6834776
[9,] 0.23512315 0.47024629 0.7648769
[10,] 0.18617917 0.37235835 0.8138208
[11,] 0.14402691 0.28805383 0.8559731
[12,] 0.16333545 0.32667090 0.8366645
[13,] 0.11476063 0.22952125 0.8852394
[14,] 0.09978389 0.19956779 0.9002161
[15,] 0.07592121 0.15184242 0.9240788
[16,] 0.09211651 0.18423301 0.9078835
[17,] 0.07569894 0.15139789 0.9243011
[18,] 0.05257400 0.10514799 0.9474260
[19,] 0.04989644 0.09979289 0.9501036
[20,] 0.06888612 0.13777224 0.9311139
[21,] 0.05525382 0.11050764 0.9447462
[22,] 0.04251682 0.08503364 0.9574832
[23,] 0.03479878 0.06959755 0.9652012
[24,] 0.02274592 0.04549183 0.9772541
[25,] 0.02419006 0.04838012 0.9758099
[26,] 0.02272012 0.04544024 0.9772799
[27,] 0.01755541 0.03511081 0.9824446
[28,] 0.01813913 0.03627827 0.9818609
[29,] 0.02659588 0.05319177 0.9734041
[30,] 0.01698376 0.03396752 0.9830162
[31,] 0.02187212 0.04374424 0.9781279
[32,] 0.07779808 0.15559616 0.9222019
[33,] 0.06082651 0.12165303 0.9391735
[34,] 0.04267438 0.08534876 0.9573256
[35,] 0.02816617 0.05633234 0.9718338
[36,] 0.04402366 0.08804732 0.9559763
[37,] 0.06533810 0.13067620 0.9346619
[38,] 0.04908917 0.09817833 0.9509108
[39,] 0.07754726 0.15509451 0.9224527
[40,] 0.06758464 0.13516928 0.9324154
[41,] 0.04745001 0.09490001 0.9525500
[42,] 0.21291202 0.42582404 0.7870880
[43,] 0.26626130 0.53252260 0.7337387
[44,] 0.20273962 0.40547924 0.7972604
[45,] 0.47826136 0.95652273 0.5217386
[46,] 0.47954822 0.95909643 0.5204518
[47,] 0.44340885 0.88681771 0.5565911
[48,] 0.55117510 0.89764981 0.4488249
[49,] 0.57601804 0.84796392 0.4239820
[50,] 0.46376576 0.92753152 0.5362342
[51,] 0.48913101 0.97826201 0.5108690
> postscript(file="/var/www/html/rcomp/tmp/1zrjh1229279594.ps",horizontal=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/2mw8b1229279594.ps",horizontal=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/3uywy1229279594.ps",horizontal=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/4dj2g1229279594.ps",horizontal=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/591dc1229279594.ps",horizontal=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 = 84
Frequency = 1
1 2 3 4 5 6
206.8095238 224.1380952 169.0809524 -20.0476190 79.6238095 -170.3190476
7 8 9 10 11 12
-36.4047619 29.1238095 -112.0904762 -63.5190476 -191.4904762 -114.9047619
13 14 15 16 17 18
64.7674603 11.1960317 -173.5611111 62.8103175 12.4817460 -18.3611111
19 20 21 22 23 24
170.8531746 -111.5182540 -60.9325397 100.5388889 -70.8325397 -106.0468254
25 26 27 28 29 30
-123.7334921 -0.7049206 -98.9620635 95.2093651 -12.1192063 -39.1620635
31 32 33 34 35 36
85.0522222 -232.6192063 -90.6334921 17.2379365 -58.4334921 99.9522222
37 38 39 40 41 42
-93.6344444 -86.2058730 23.0369841 30.1084127 -150.8201587 2.0369841
43 44 45 46 47 48
106.1512698 6.8798413 155.3655556 44.9369841 152.8655556 325.2512698
49 50 51 52 53 54
-50.6353968 6.6931746 15.2360317 186.7074603 -153.2211111 24.2360317
55 56 57 58 59 60
-200.4496825 47.2788889 38.5646032 -236.6639683 8.2646032 32.2503175
61 62 63 64 65 66
-125.9363492 -73.0077778 63.3350794 -177.0934921 152.3779365 236.0350794
67 68 69 70 71 72
-18.2506349 254.2779365 215.9636508 122.4350794 265.9636508 128.1493651
73 74 75 76 77 78
122.3626984 -82.1087302 1.8341270 -177.6944444 71.6769841 -34.4658730
79 80 81 82 83 84
-106.9515873 6.5769841 -146.2373016 15.0341270 -106.3373016 -364.6515873
> postscript(file="/var/www/html/rcomp/tmp/6p3vy1229279594.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 84
Frequency = 1
lag(myerror, k = 1) myerror
0 206.8095238 NA
1 224.1380952 206.8095238
2 169.0809524 224.1380952
3 -20.0476190 169.0809524
4 79.6238095 -20.0476190
5 -170.3190476 79.6238095
6 -36.4047619 -170.3190476
7 29.1238095 -36.4047619
8 -112.0904762 29.1238095
9 -63.5190476 -112.0904762
10 -191.4904762 -63.5190476
11 -114.9047619 -191.4904762
12 64.7674603 -114.9047619
13 11.1960317 64.7674603
14 -173.5611111 11.1960317
15 62.8103175 -173.5611111
16 12.4817460 62.8103175
17 -18.3611111 12.4817460
18 170.8531746 -18.3611111
19 -111.5182540 170.8531746
20 -60.9325397 -111.5182540
21 100.5388889 -60.9325397
22 -70.8325397 100.5388889
23 -106.0468254 -70.8325397
24 -123.7334921 -106.0468254
25 -0.7049206 -123.7334921
26 -98.9620635 -0.7049206
27 95.2093651 -98.9620635
28 -12.1192063 95.2093651
29 -39.1620635 -12.1192063
30 85.0522222 -39.1620635
31 -232.6192063 85.0522222
32 -90.6334921 -232.6192063
33 17.2379365 -90.6334921
34 -58.4334921 17.2379365
35 99.9522222 -58.4334921
36 -93.6344444 99.9522222
37 -86.2058730 -93.6344444
38 23.0369841 -86.2058730
39 30.1084127 23.0369841
40 -150.8201587 30.1084127
41 2.0369841 -150.8201587
42 106.1512698 2.0369841
43 6.8798413 106.1512698
44 155.3655556 6.8798413
45 44.9369841 155.3655556
46 152.8655556 44.9369841
47 325.2512698 152.8655556
48 -50.6353968 325.2512698
49 6.6931746 -50.6353968
50 15.2360317 6.6931746
51 186.7074603 15.2360317
52 -153.2211111 186.7074603
53 24.2360317 -153.2211111
54 -200.4496825 24.2360317
55 47.2788889 -200.4496825
56 38.5646032 47.2788889
57 -236.6639683 38.5646032
58 8.2646032 -236.6639683
59 32.2503175 8.2646032
60 -125.9363492 32.2503175
61 -73.0077778 -125.9363492
62 63.3350794 -73.0077778
63 -177.0934921 63.3350794
64 152.3779365 -177.0934921
65 236.0350794 152.3779365
66 -18.2506349 236.0350794
67 254.2779365 -18.2506349
68 215.9636508 254.2779365
69 122.4350794 215.9636508
70 265.9636508 122.4350794
71 128.1493651 265.9636508
72 122.3626984 128.1493651
73 -82.1087302 122.3626984
74 1.8341270 -82.1087302
75 -177.6944444 1.8341270
76 71.6769841 -177.6944444
77 -34.4658730 71.6769841
78 -106.9515873 -34.4658730
79 6.5769841 -106.9515873
80 -146.2373016 6.5769841
81 15.0341270 -146.2373016
82 -106.3373016 15.0341270
83 -364.6515873 -106.3373016
84 NA -364.6515873
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 224.1380952 206.8095238
[2,] 169.0809524 224.1380952
[3,] -20.0476190 169.0809524
[4,] 79.6238095 -20.0476190
[5,] -170.3190476 79.6238095
[6,] -36.4047619 -170.3190476
[7,] 29.1238095 -36.4047619
[8,] -112.0904762 29.1238095
[9,] -63.5190476 -112.0904762
[10,] -191.4904762 -63.5190476
[11,] -114.9047619 -191.4904762
[12,] 64.7674603 -114.9047619
[13,] 11.1960317 64.7674603
[14,] -173.5611111 11.1960317
[15,] 62.8103175 -173.5611111
[16,] 12.4817460 62.8103175
[17,] -18.3611111 12.4817460
[18,] 170.8531746 -18.3611111
[19,] -111.5182540 170.8531746
[20,] -60.9325397 -111.5182540
[21,] 100.5388889 -60.9325397
[22,] -70.8325397 100.5388889
[23,] -106.0468254 -70.8325397
[24,] -123.7334921 -106.0468254
[25,] -0.7049206 -123.7334921
[26,] -98.9620635 -0.7049206
[27,] 95.2093651 -98.9620635
[28,] -12.1192063 95.2093651
[29,] -39.1620635 -12.1192063
[30,] 85.0522222 -39.1620635
[31,] -232.6192063 85.0522222
[32,] -90.6334921 -232.6192063
[33,] 17.2379365 -90.6334921
[34,] -58.4334921 17.2379365
[35,] 99.9522222 -58.4334921
[36,] -93.6344444 99.9522222
[37,] -86.2058730 -93.6344444
[38,] 23.0369841 -86.2058730
[39,] 30.1084127 23.0369841
[40,] -150.8201587 30.1084127
[41,] 2.0369841 -150.8201587
[42,] 106.1512698 2.0369841
[43,] 6.8798413 106.1512698
[44,] 155.3655556 6.8798413
[45,] 44.9369841 155.3655556
[46,] 152.8655556 44.9369841
[47,] 325.2512698 152.8655556
[48,] -50.6353968 325.2512698
[49,] 6.6931746 -50.6353968
[50,] 15.2360317 6.6931746
[51,] 186.7074603 15.2360317
[52,] -153.2211111 186.7074603
[53,] 24.2360317 -153.2211111
[54,] -200.4496825 24.2360317
[55,] 47.2788889 -200.4496825
[56,] 38.5646032 47.2788889
[57,] -236.6639683 38.5646032
[58,] 8.2646032 -236.6639683
[59,] 32.2503175 8.2646032
[60,] -125.9363492 32.2503175
[61,] -73.0077778 -125.9363492
[62,] 63.3350794 -73.0077778
[63,] -177.0934921 63.3350794
[64,] 152.3779365 -177.0934921
[65,] 236.0350794 152.3779365
[66,] -18.2506349 236.0350794
[67,] 254.2779365 -18.2506349
[68,] 215.9636508 254.2779365
[69,] 122.4350794 215.9636508
[70,] 265.9636508 122.4350794
[71,] 128.1493651 265.9636508
[72,] 122.3626984 128.1493651
[73,] -82.1087302 122.3626984
[74,] 1.8341270 -82.1087302
[75,] -177.6944444 1.8341270
[76,] 71.6769841 -177.6944444
[77,] -34.4658730 71.6769841
[78,] -106.9515873 -34.4658730
[79,] 6.5769841 -106.9515873
[80,] -146.2373016 6.5769841
[81,] 15.0341270 -146.2373016
[82,] -106.3373016 15.0341270
[83,] -364.6515873 -106.3373016
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 224.1380952 206.8095238
2 169.0809524 224.1380952
3 -20.0476190 169.0809524
4 79.6238095 -20.0476190
5 -170.3190476 79.6238095
6 -36.4047619 -170.3190476
7 29.1238095 -36.4047619
8 -112.0904762 29.1238095
9 -63.5190476 -112.0904762
10 -191.4904762 -63.5190476
11 -114.9047619 -191.4904762
12 64.7674603 -114.9047619
13 11.1960317 64.7674603
14 -173.5611111 11.1960317
15 62.8103175 -173.5611111
16 12.4817460 62.8103175
17 -18.3611111 12.4817460
18 170.8531746 -18.3611111
19 -111.5182540 170.8531746
20 -60.9325397 -111.5182540
21 100.5388889 -60.9325397
22 -70.8325397 100.5388889
23 -106.0468254 -70.8325397
24 -123.7334921 -106.0468254
25 -0.7049206 -123.7334921
26 -98.9620635 -0.7049206
27 95.2093651 -98.9620635
28 -12.1192063 95.2093651
29 -39.1620635 -12.1192063
30 85.0522222 -39.1620635
31 -232.6192063 85.0522222
32 -90.6334921 -232.6192063
33 17.2379365 -90.6334921
34 -58.4334921 17.2379365
35 99.9522222 -58.4334921
36 -93.6344444 99.9522222
37 -86.2058730 -93.6344444
38 23.0369841 -86.2058730
39 30.1084127 23.0369841
40 -150.8201587 30.1084127
41 2.0369841 -150.8201587
42 106.1512698 2.0369841
43 6.8798413 106.1512698
44 155.3655556 6.8798413
45 44.9369841 155.3655556
46 152.8655556 44.9369841
47 325.2512698 152.8655556
48 -50.6353968 325.2512698
49 6.6931746 -50.6353968
50 15.2360317 6.6931746
51 186.7074603 15.2360317
52 -153.2211111 186.7074603
53 24.2360317 -153.2211111
54 -200.4496825 24.2360317
55 47.2788889 -200.4496825
56 38.5646032 47.2788889
57 -236.6639683 38.5646032
58 8.2646032 -236.6639683
59 32.2503175 8.2646032
60 -125.9363492 32.2503175
61 -73.0077778 -125.9363492
62 63.3350794 -73.0077778
63 -177.0934921 63.3350794
64 152.3779365 -177.0934921
65 236.0350794 152.3779365
66 -18.2506349 236.0350794
67 254.2779365 -18.2506349
68 215.9636508 254.2779365
69 122.4350794 215.9636508
70 265.9636508 122.4350794
71 128.1493651 265.9636508
72 122.3626984 128.1493651
73 -82.1087302 122.3626984
74 1.8341270 -82.1087302
75 -177.6944444 1.8341270
76 71.6769841 -177.6944444
77 -34.4658730 71.6769841
78 -106.9515873 -34.4658730
79 6.5769841 -106.9515873
80 -146.2373016 6.5769841
81 15.0341270 -146.2373016
82 -106.3373016 15.0341270
83 -364.6515873 -106.3373016
> 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/78e311229279594.ps",horizontal=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/8g3vt1229279594.ps",horizontal=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/9i8381229279594.ps",horizontal=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/10rock1229279594.ps",horizontal=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/1174lt1229279594.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/12ii1r1229279594.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/13imnv1229279594.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/14k9dn1229279594.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/15oqs71229279594.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/16n5st1229279595.tab")
+ }
>
> system("convert tmp/1zrjh1229279594.ps tmp/1zrjh1229279594.png")
> system("convert tmp/2mw8b1229279594.ps tmp/2mw8b1229279594.png")
> system("convert tmp/3uywy1229279594.ps tmp/3uywy1229279594.png")
> system("convert tmp/4dj2g1229279594.ps tmp/4dj2g1229279594.png")
> system("convert tmp/591dc1229279594.ps tmp/591dc1229279594.png")
> system("convert tmp/6p3vy1229279594.ps tmp/6p3vy1229279594.png")
> system("convert tmp/78e311229279594.ps tmp/78e311229279594.png")
> system("convert tmp/8g3vt1229279594.ps tmp/8g3vt1229279594.png")
> system("convert tmp/9i8381229279594.ps tmp/9i8381229279594.png")
> system("convert tmp/10rock1229279594.ps tmp/10rock1229279594.png")
>
>
> proc.time()
user system elapsed
7.184 2.805 7.621