R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-pc-linux-gnu (32-bit)
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(111
+ ,45
+ ,27
+ ,52
+ ,102
+ ,50
+ ,18
+ ,55
+ ,108
+ ,49
+ ,19
+ ,80
+ ,109
+ ,55
+ ,20
+ ,45
+ ,118
+ ,39
+ ,29
+ ,60
+ ,79
+ ,68
+ ,46
+ ,34
+ ,88
+ ,69
+ ,27
+ ,45
+ ,102
+ ,56
+ ,23
+ ,68
+ ,105
+ ,58
+ ,29
+ ,26
+ ,92
+ ,48
+ ,38
+ ,70
+ ,131
+ ,34
+ ,20
+ ,85
+ ,104
+ ,50
+ ,37
+ ,54
+ ,83
+ ,76
+ ,32
+ ,55
+ ,84
+ ,49
+ ,26
+ ,40
+ ,85
+ ,51
+ ,40
+ ,55
+ ,110
+ ,53
+ ,30
+ ,50
+ ,121
+ ,36
+ ,26
+ ,71
+ ,120
+ ,62
+ ,23
+ ,55
+ ,100
+ ,46
+ ,27
+ ,70
+ ,94
+ ,50
+ ,38
+ ,55
+ ,89
+ ,47
+ ,25
+ ,60
+ ,93
+ ,50
+ ,33
+ ,65
+ ,128
+ ,44
+ ,45
+ ,66
+ ,84
+ ,50
+ ,34
+ ,55
+ ,127
+ ,29
+ ,20
+ ,90
+ ,106
+ ,49
+ ,24
+ ,55
+ ,129
+ ,26
+ ,26
+ ,60
+ ,82
+ ,79
+ ,26
+ ,35
+ ,106
+ ,53
+ ,39
+ ,55
+ ,109
+ ,53
+ ,27
+ ,26
+ ,91
+ ,72
+ ,18
+ ,14
+ ,111
+ ,35
+ ,34
+ ,45
+ ,105
+ ,42
+ ,25
+ ,35
+ ,118
+ ,37
+ ,26
+ ,65
+ ,103
+ ,46
+ ,28
+ ,35
+ ,101
+ ,48
+ ,21
+ ,60
+ ,101
+ ,46
+ ,39
+ ,60
+ ,95
+ ,49
+ ,25
+ ,60
+ ,108
+ ,65
+ ,29
+ ,65
+ ,95
+ ,52
+ ,37
+ ,45
+ ,98
+ ,75
+ ,34
+ ,20
+ ,82
+ ,58
+ ,30
+ ,50
+ ,100
+ ,43
+ ,28
+ ,60
+ ,100
+ ,60
+ ,25
+ ,48
+ ,107
+ ,43
+ ,27
+ ,40
+ ,95
+ ,51
+ ,33
+ ,55
+ ,97
+ ,70
+ ,30
+ ,54
+ ,93
+ ,69
+ ,26
+ ,40
+ ,81
+ ,65
+ ,18
+ ,40
+ ,89
+ ,63
+ ,21
+ ,34
+ ,111
+ ,44
+ ,39
+ ,60
+ ,95
+ ,61
+ ,36
+ ,30
+ ,106
+ ,40
+ ,32
+ ,75
+ ,83
+ ,62
+ ,23
+ ,24
+ ,81
+ ,59
+ ,27
+ ,30
+ ,115
+ ,47
+ ,45
+ ,80
+ ,112
+ ,50
+ ,24
+ ,60
+ ,92
+ ,50
+ ,29
+ ,46
+ ,85
+ ,65
+ ,21
+ ,35
+ ,95
+ ,54
+ ,28
+ ,60
+ ,115
+ ,44
+ ,37
+ ,75
+ ,91
+ ,66
+ ,22
+ ,54
+ ,107
+ ,34
+ ,31
+ ,78
+ ,102
+ ,74
+ ,32
+ ,20
+ ,86
+ ,57
+ ,20
+ ,45
+ ,96
+ ,60
+ ,33
+ ,60
+ ,114
+ ,36
+ ,32
+ ,70
+ ,105
+ ,50
+ ,18
+ ,35
+ ,82
+ ,60
+ ,44
+ ,20
+ ,120
+ ,45
+ ,24
+ ,60
+ ,88
+ ,55
+ ,21
+ ,20
+ ,90
+ ,44
+ ,29
+ ,50
+ ,85
+ ,57
+ ,30
+ ,50
+ ,106
+ ,33
+ ,37
+ ,75
+ ,109
+ ,30
+ ,33
+ ,70
+ ,75
+ ,64
+ ,25
+ ,20
+ ,91
+ ,49
+ ,19
+ ,45
+ ,96
+ ,76
+ ,16
+ ,20
+ ,108
+ ,40
+ ,31
+ ,50
+ ,86
+ ,48
+ ,29
+ ,55
+ ,98
+ ,65
+ ,37
+ ,15
+ ,99
+ ,50
+ ,41
+ ,26
+ ,95
+ ,70
+ ,28
+ ,25
+ ,88
+ ,78
+ ,19
+ ,30
+ ,111
+ ,44
+ ,28
+ ,60
+ ,103
+ ,48
+ ,33
+ ,40
+ ,107
+ ,52
+ ,32
+ ,40
+ ,118
+ ,40
+ ,28
+ ,50)
+ ,dim=c(4
+ ,88)
+ ,dimnames=list(c('IQ'
+ ,'Add'
+ ,'MumAge'
+ ,'Grade')
+ ,1:88))
> y <- array(NA,dim=c(4,88),dimnames=list(c('IQ','Add','MumAge','Grade'),1:88))
> 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 = 'Do not include Seasonal Dummies'
> par1 = '1'
> par3 <- 'No Linear Trend'
> par2 <- 'Do not include Seasonal 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, 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
IQ Add MumAge Grade
1 111 45 27 52
2 102 50 18 55
3 108 49 19 80
4 109 55 20 45
5 118 39 29 60
6 79 68 46 34
7 88 69 27 45
8 102 56 23 68
9 105 58 29 26
10 92 48 38 70
11 131 34 20 85
12 104 50 37 54
13 83 76 32 55
14 84 49 26 40
15 85 51 40 55
16 110 53 30 50
17 121 36 26 71
18 120 62 23 55
19 100 46 27 70
20 94 50 38 55
21 89 47 25 60
22 93 50 33 65
23 128 44 45 66
24 84 50 34 55
25 127 29 20 90
26 106 49 24 55
27 129 26 26 60
28 82 79 26 35
29 106 53 39 55
30 109 53 27 26
31 91 72 18 14
32 111 35 34 45
33 105 42 25 35
34 118 37 26 65
35 103 46 28 35
36 101 48 21 60
37 101 46 39 60
38 95 49 25 60
39 108 65 29 65
40 95 52 37 45
41 98 75 34 20
42 82 58 30 50
43 100 43 28 60
44 100 60 25 48
45 107 43 27 40
46 95 51 33 55
47 97 70 30 54
48 93 69 26 40
49 81 65 18 40
50 89 63 21 34
51 111 44 39 60
52 95 61 36 30
53 106 40 32 75
54 83 62 23 24
55 81 59 27 30
56 115 47 45 80
57 112 50 24 60
58 92 50 29 46
59 85 65 21 35
60 95 54 28 60
61 115 44 37 75
62 91 66 22 54
63 107 34 31 78
64 102 74 32 20
65 86 57 20 45
66 96 60 33 60
67 114 36 32 70
68 105 50 18 35
69 82 60 44 20
70 120 45 24 60
71 88 55 21 20
72 90 44 29 50
73 85 57 30 50
74 106 33 37 75
75 109 30 33 70
76 75 64 25 20
77 91 49 19 45
78 96 76 16 20
79 108 40 31 50
80 86 48 29 55
81 98 65 37 15
82 99 50 41 26
83 95 70 28 25
84 88 78 19 30
85 111 44 28 60
86 103 48 33 40
87 107 52 32 40
88 118 40 28 50
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Add MumAge Grade
124.8063 -0.5396 -0.1256 0.1465
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-17.613 -7.252 1.109 7.278 23.480
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 124.80627 10.32309 12.090 < 2e-16 ***
Add -0.53956 0.11539 -4.676 1.11e-05 ***
MumAge -0.12565 0.15235 -0.825 0.4119
Grade 0.14648 0.07884 1.858 0.0667 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.852 on 84 degrees of freedom
Multiple R-squared: 0.4236, Adjusted R-squared: 0.403
F-statistic: 20.58 on 3 and 84 DF, p-value: 4.349e-10
> 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.16350387 0.32700774 0.83649613
[2,] 0.07381400 0.14762799 0.92618600
[3,] 0.04602220 0.09204439 0.95397780
[4,] 0.03808568 0.07617135 0.96191432
[5,] 0.06618619 0.13237238 0.93381381
[6,] 0.03940005 0.07880009 0.96059995
[7,] 0.03610639 0.07221278 0.96389361
[8,] 0.41279139 0.82558278 0.58720861
[9,] 0.44968857 0.89937713 0.55031143
[10,] 0.46412823 0.92825647 0.53587177
[11,] 0.40378970 0.80757941 0.59621030
[12,] 0.71345958 0.57308085 0.28654042
[13,] 0.69072365 0.61855270 0.30927635
[14,] 0.63252072 0.73495857 0.36747928
[15,] 0.78666127 0.42667746 0.21333873
[16,] 0.76514910 0.46970181 0.23485090
[17,] 0.96356754 0.07286491 0.03643246
[18,] 0.98226592 0.03546817 0.01773408
[19,] 0.97804468 0.04391064 0.02195532
[20,] 0.96824547 0.06350907 0.03175453
[21,] 0.97293584 0.05412832 0.02706416
[22,] 0.96212813 0.07574375 0.03787187
[23,] 0.95532135 0.08935730 0.04467865
[24,] 0.95887203 0.08225594 0.04112797
[25,] 0.94582834 0.10834331 0.05417166
[26,] 0.92914797 0.14170406 0.07085203
[27,] 0.91361390 0.17277221 0.08638610
[28,] 0.90750137 0.18499726 0.09249863
[29,] 0.88419089 0.23161823 0.11580911
[30,] 0.86208902 0.27582196 0.13791098
[31,] 0.82790548 0.34418905 0.17209452
[32,] 0.82162960 0.35674080 0.17837040
[33,] 0.84663354 0.30673292 0.15336646
[34,] 0.81400607 0.37198786 0.18599393
[35,] 0.85004880 0.29990240 0.14995120
[36,] 0.90021883 0.19956235 0.09978117
[37,] 0.88334042 0.23331917 0.11665958
[38,] 0.85445583 0.29108834 0.14554417
[39,] 0.82909349 0.34181301 0.17090651
[40,] 0.80482402 0.39035197 0.19517598
[41,] 0.76889840 0.46220319 0.23110160
[42,] 0.71878948 0.56242103 0.28121052
[43,] 0.75028124 0.49943752 0.24971876
[44,] 0.70641064 0.58717873 0.29358936
[45,] 0.66834595 0.66330810 0.33165405
[46,] 0.61053762 0.77892476 0.38946238
[47,] 0.55508943 0.88982115 0.44491057
[48,] 0.53862423 0.92275153 0.46137577
[49,] 0.58127028 0.83745945 0.41872972
[50,] 0.55775214 0.88449572 0.44224786
[51,] 0.55678754 0.88642492 0.44321246
[52,] 0.53497800 0.93004400 0.46502200
[53,] 0.50024150 0.99951701 0.49975850
[54,] 0.45288086 0.90576173 0.54711914
[55,] 0.43377846 0.86755691 0.56622154
[56,] 0.37249979 0.74499958 0.62750021
[57,] 0.32498369 0.64996737 0.67501631
[58,] 0.44007440 0.88014881 0.55992560
[59,] 0.46534273 0.93068546 0.53465727
[60,] 0.39134192 0.78268384 0.60865808
[61,] 0.33579448 0.67158896 0.66420552
[62,] 0.28211003 0.56422005 0.71788997
[63,] 0.27264417 0.54528834 0.72735583
[64,] 0.42415828 0.84831657 0.57584172
[65,] 0.38340709 0.76681417 0.61659291
[66,] 0.43158292 0.86316583 0.56841708
[67,] 0.46403805 0.92807610 0.53596195
[68,] 0.39872881 0.79745762 0.60127119
[69,] 0.31745694 0.63491389 0.68254306
[70,] 0.58251636 0.83496727 0.41748364
[71,] 0.80503481 0.38993039 0.19496519
[72,] 0.77279667 0.45440665 0.22720333
[73,] 0.68189895 0.63620210 0.31810105
[74,] 0.94711225 0.10577549 0.05288775
[75,] 0.93028107 0.13943786 0.06971893
> postscript(file="/var/wessaorg/rcomp/tmp/1lzml1356040516.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/wessaorg/rcomp/tmp/2po151356040516.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/wessaorg/rcomp/tmp/3k4eq1356040516.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/wessaorg/rcomp/tmp/4khvy1356040516.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/wessaorg/rcomp/tmp/5kh9w1356040516.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 = 88
Frequency = 1
1 2 3 4 5 6
6.2495150 -1.6229716 0.3012408 9.7908542 9.0916787 -8.3168367
7 8 9 10 11 12
-2.7758321 0.3384248 11.3233867 -12.3862391 14.6011843 2.9108267
13 14 15 16 17 18
-4.8354547 -16.9602097 -15.3191476 10.2358533 8.4848396 23.4799344
19 20 21 22 23 24
-6.8474840 -7.1100000 -16.0944735 -10.2029953 22.9209809 -17.6125943
25 26 27 28 29 30
7.1710318 2.5913646 12.7005146 -2.0411757 6.6343144 12.3743133
31 32 33 34 35 36
5.2527288 2.7588292 0.8696311 6.9032463 1.4047979 -4.0575126
37 38 39 40 41 42
-2.8749486 -9.0153630 12.3877393 -3.6917857 15.0029205 -15.0663704
43 44 45 46 47 48
-6.8757488 3.6774477 2.9281073 -6.1986876 5.8223918 2.8308955
49 50 51 52 53 54
-12.3325142 -4.1558275 6.0459408 3.2356916 -4.1889488 -8.9793333
55 56 57 58 59 60
-12.9742561 9.4889934 8.3985437 -8.9225601 -7.2231922 -5.9406410
61 62 63 64 65 66
7.5975151 -3.3410179 -6.9913547 18.2120680 -12.1300353 -1.0750665
67 68 69 70 71 72
2.3852063 4.3065331 -7.8339226 13.7007674 -7.4216163 -14.7457926
73 74 75 76 77 78
-12.6059256 -7.3375927 -5.7264767 -15.0630246 -11.5721259 11.2808013
79 80 81 82 83 84
1.3472835 -17.3199478 10.7166898 2.5147276 7.8188765 3.2721052
85 86 87 88
4.6638064 2.3797751 8.4123476 10.9703378
> postscript(file="/var/wessaorg/rcomp/tmp/6opa91356040516.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 = 88
Frequency = 1
lag(myerror, k = 1) myerror
0 6.2495150 NA
1 -1.6229716 6.2495150
2 0.3012408 -1.6229716
3 9.7908542 0.3012408
4 9.0916787 9.7908542
5 -8.3168367 9.0916787
6 -2.7758321 -8.3168367
7 0.3384248 -2.7758321
8 11.3233867 0.3384248
9 -12.3862391 11.3233867
10 14.6011843 -12.3862391
11 2.9108267 14.6011843
12 -4.8354547 2.9108267
13 -16.9602097 -4.8354547
14 -15.3191476 -16.9602097
15 10.2358533 -15.3191476
16 8.4848396 10.2358533
17 23.4799344 8.4848396
18 -6.8474840 23.4799344
19 -7.1100000 -6.8474840
20 -16.0944735 -7.1100000
21 -10.2029953 -16.0944735
22 22.9209809 -10.2029953
23 -17.6125943 22.9209809
24 7.1710318 -17.6125943
25 2.5913646 7.1710318
26 12.7005146 2.5913646
27 -2.0411757 12.7005146
28 6.6343144 -2.0411757
29 12.3743133 6.6343144
30 5.2527288 12.3743133
31 2.7588292 5.2527288
32 0.8696311 2.7588292
33 6.9032463 0.8696311
34 1.4047979 6.9032463
35 -4.0575126 1.4047979
36 -2.8749486 -4.0575126
37 -9.0153630 -2.8749486
38 12.3877393 -9.0153630
39 -3.6917857 12.3877393
40 15.0029205 -3.6917857
41 -15.0663704 15.0029205
42 -6.8757488 -15.0663704
43 3.6774477 -6.8757488
44 2.9281073 3.6774477
45 -6.1986876 2.9281073
46 5.8223918 -6.1986876
47 2.8308955 5.8223918
48 -12.3325142 2.8308955
49 -4.1558275 -12.3325142
50 6.0459408 -4.1558275
51 3.2356916 6.0459408
52 -4.1889488 3.2356916
53 -8.9793333 -4.1889488
54 -12.9742561 -8.9793333
55 9.4889934 -12.9742561
56 8.3985437 9.4889934
57 -8.9225601 8.3985437
58 -7.2231922 -8.9225601
59 -5.9406410 -7.2231922
60 7.5975151 -5.9406410
61 -3.3410179 7.5975151
62 -6.9913547 -3.3410179
63 18.2120680 -6.9913547
64 -12.1300353 18.2120680
65 -1.0750665 -12.1300353
66 2.3852063 -1.0750665
67 4.3065331 2.3852063
68 -7.8339226 4.3065331
69 13.7007674 -7.8339226
70 -7.4216163 13.7007674
71 -14.7457926 -7.4216163
72 -12.6059256 -14.7457926
73 -7.3375927 -12.6059256
74 -5.7264767 -7.3375927
75 -15.0630246 -5.7264767
76 -11.5721259 -15.0630246
77 11.2808013 -11.5721259
78 1.3472835 11.2808013
79 -17.3199478 1.3472835
80 10.7166898 -17.3199478
81 2.5147276 10.7166898
82 7.8188765 2.5147276
83 3.2721052 7.8188765
84 4.6638064 3.2721052
85 2.3797751 4.6638064
86 8.4123476 2.3797751
87 10.9703378 8.4123476
88 NA 10.9703378
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.6229716 6.2495150
[2,] 0.3012408 -1.6229716
[3,] 9.7908542 0.3012408
[4,] 9.0916787 9.7908542
[5,] -8.3168367 9.0916787
[6,] -2.7758321 -8.3168367
[7,] 0.3384248 -2.7758321
[8,] 11.3233867 0.3384248
[9,] -12.3862391 11.3233867
[10,] 14.6011843 -12.3862391
[11,] 2.9108267 14.6011843
[12,] -4.8354547 2.9108267
[13,] -16.9602097 -4.8354547
[14,] -15.3191476 -16.9602097
[15,] 10.2358533 -15.3191476
[16,] 8.4848396 10.2358533
[17,] 23.4799344 8.4848396
[18,] -6.8474840 23.4799344
[19,] -7.1100000 -6.8474840
[20,] -16.0944735 -7.1100000
[21,] -10.2029953 -16.0944735
[22,] 22.9209809 -10.2029953
[23,] -17.6125943 22.9209809
[24,] 7.1710318 -17.6125943
[25,] 2.5913646 7.1710318
[26,] 12.7005146 2.5913646
[27,] -2.0411757 12.7005146
[28,] 6.6343144 -2.0411757
[29,] 12.3743133 6.6343144
[30,] 5.2527288 12.3743133
[31,] 2.7588292 5.2527288
[32,] 0.8696311 2.7588292
[33,] 6.9032463 0.8696311
[34,] 1.4047979 6.9032463
[35,] -4.0575126 1.4047979
[36,] -2.8749486 -4.0575126
[37,] -9.0153630 -2.8749486
[38,] 12.3877393 -9.0153630
[39,] -3.6917857 12.3877393
[40,] 15.0029205 -3.6917857
[41,] -15.0663704 15.0029205
[42,] -6.8757488 -15.0663704
[43,] 3.6774477 -6.8757488
[44,] 2.9281073 3.6774477
[45,] -6.1986876 2.9281073
[46,] 5.8223918 -6.1986876
[47,] 2.8308955 5.8223918
[48,] -12.3325142 2.8308955
[49,] -4.1558275 -12.3325142
[50,] 6.0459408 -4.1558275
[51,] 3.2356916 6.0459408
[52,] -4.1889488 3.2356916
[53,] -8.9793333 -4.1889488
[54,] -12.9742561 -8.9793333
[55,] 9.4889934 -12.9742561
[56,] 8.3985437 9.4889934
[57,] -8.9225601 8.3985437
[58,] -7.2231922 -8.9225601
[59,] -5.9406410 -7.2231922
[60,] 7.5975151 -5.9406410
[61,] -3.3410179 7.5975151
[62,] -6.9913547 -3.3410179
[63,] 18.2120680 -6.9913547
[64,] -12.1300353 18.2120680
[65,] -1.0750665 -12.1300353
[66,] 2.3852063 -1.0750665
[67,] 4.3065331 2.3852063
[68,] -7.8339226 4.3065331
[69,] 13.7007674 -7.8339226
[70,] -7.4216163 13.7007674
[71,] -14.7457926 -7.4216163
[72,] -12.6059256 -14.7457926
[73,] -7.3375927 -12.6059256
[74,] -5.7264767 -7.3375927
[75,] -15.0630246 -5.7264767
[76,] -11.5721259 -15.0630246
[77,] 11.2808013 -11.5721259
[78,] 1.3472835 11.2808013
[79,] -17.3199478 1.3472835
[80,] 10.7166898 -17.3199478
[81,] 2.5147276 10.7166898
[82,] 7.8188765 2.5147276
[83,] 3.2721052 7.8188765
[84,] 4.6638064 3.2721052
[85,] 2.3797751 4.6638064
[86,] 8.4123476 2.3797751
[87,] 10.9703378 8.4123476
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.6229716 6.2495150
2 0.3012408 -1.6229716
3 9.7908542 0.3012408
4 9.0916787 9.7908542
5 -8.3168367 9.0916787
6 -2.7758321 -8.3168367
7 0.3384248 -2.7758321
8 11.3233867 0.3384248
9 -12.3862391 11.3233867
10 14.6011843 -12.3862391
11 2.9108267 14.6011843
12 -4.8354547 2.9108267
13 -16.9602097 -4.8354547
14 -15.3191476 -16.9602097
15 10.2358533 -15.3191476
16 8.4848396 10.2358533
17 23.4799344 8.4848396
18 -6.8474840 23.4799344
19 -7.1100000 -6.8474840
20 -16.0944735 -7.1100000
21 -10.2029953 -16.0944735
22 22.9209809 -10.2029953
23 -17.6125943 22.9209809
24 7.1710318 -17.6125943
25 2.5913646 7.1710318
26 12.7005146 2.5913646
27 -2.0411757 12.7005146
28 6.6343144 -2.0411757
29 12.3743133 6.6343144
30 5.2527288 12.3743133
31 2.7588292 5.2527288
32 0.8696311 2.7588292
33 6.9032463 0.8696311
34 1.4047979 6.9032463
35 -4.0575126 1.4047979
36 -2.8749486 -4.0575126
37 -9.0153630 -2.8749486
38 12.3877393 -9.0153630
39 -3.6917857 12.3877393
40 15.0029205 -3.6917857
41 -15.0663704 15.0029205
42 -6.8757488 -15.0663704
43 3.6774477 -6.8757488
44 2.9281073 3.6774477
45 -6.1986876 2.9281073
46 5.8223918 -6.1986876
47 2.8308955 5.8223918
48 -12.3325142 2.8308955
49 -4.1558275 -12.3325142
50 6.0459408 -4.1558275
51 3.2356916 6.0459408
52 -4.1889488 3.2356916
53 -8.9793333 -4.1889488
54 -12.9742561 -8.9793333
55 9.4889934 -12.9742561
56 8.3985437 9.4889934
57 -8.9225601 8.3985437
58 -7.2231922 -8.9225601
59 -5.9406410 -7.2231922
60 7.5975151 -5.9406410
61 -3.3410179 7.5975151
62 -6.9913547 -3.3410179
63 18.2120680 -6.9913547
64 -12.1300353 18.2120680
65 -1.0750665 -12.1300353
66 2.3852063 -1.0750665
67 4.3065331 2.3852063
68 -7.8339226 4.3065331
69 13.7007674 -7.8339226
70 -7.4216163 13.7007674
71 -14.7457926 -7.4216163
72 -12.6059256 -14.7457926
73 -7.3375927 -12.6059256
74 -5.7264767 -7.3375927
75 -15.0630246 -5.7264767
76 -11.5721259 -15.0630246
77 11.2808013 -11.5721259
78 1.3472835 11.2808013
79 -17.3199478 1.3472835
80 10.7166898 -17.3199478
81 2.5147276 10.7166898
82 7.8188765 2.5147276
83 3.2721052 7.8188765
84 4.6638064 3.2721052
85 2.3797751 4.6638064
86 8.4123476 2.3797751
87 10.9703378 8.4123476
> 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/wessaorg/rcomp/tmp/7uofn1356040516.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/wessaorg/rcomp/tmp/8sxnw1356040516.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/wessaorg/rcomp/tmp/9qgdh1356040516.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/wessaorg/rcomp/tmp/10ff4k1356040516.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/114bks1356040516.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/wessaorg/rcomp/tmp/12xnbd1356040517.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/wessaorg/rcomp/tmp/13brrn1356040517.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/wessaorg/rcomp/tmp/14amgr1356040517.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/wessaorg/rcomp/tmp/15g8vq1356040517.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/wessaorg/rcomp/tmp/16oqsb1356040517.tab")
+ }
>
> try(system("convert tmp/1lzml1356040516.ps tmp/1lzml1356040516.png",intern=TRUE))
character(0)
> try(system("convert tmp/2po151356040516.ps tmp/2po151356040516.png",intern=TRUE))
character(0)
> try(system("convert tmp/3k4eq1356040516.ps tmp/3k4eq1356040516.png",intern=TRUE))
character(0)
> try(system("convert tmp/4khvy1356040516.ps tmp/4khvy1356040516.png",intern=TRUE))
character(0)
> try(system("convert tmp/5kh9w1356040516.ps tmp/5kh9w1356040516.png",intern=TRUE))
character(0)
> try(system("convert tmp/6opa91356040516.ps tmp/6opa91356040516.png",intern=TRUE))
character(0)
> try(system("convert tmp/7uofn1356040516.ps tmp/7uofn1356040516.png",intern=TRUE))
character(0)
> try(system("convert tmp/8sxnw1356040516.ps tmp/8sxnw1356040516.png",intern=TRUE))
character(0)
> try(system("convert tmp/9qgdh1356040516.ps tmp/9qgdh1356040516.png",intern=TRUE))
character(0)
> try(system("convert tmp/10ff4k1356040516.ps tmp/10ff4k1356040516.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.572 1.163 7.781