R version 2.8.0 (2008-10-20)
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.
Natural language support but running in an English locale
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(9
+ ,15
+ ,6
+ ,25
+ ,68
+ ,14
+ ,10
+ ,8
+ ,23
+ ,48
+ ,8
+ ,10
+ ,7
+ ,17
+ ,44
+ ,8
+ ,12
+ ,9
+ ,19
+ ,67
+ ,14
+ ,9
+ ,8
+ ,29
+ ,46
+ ,15
+ ,18
+ ,11
+ ,23
+ ,54
+ ,9
+ ,14
+ ,9
+ ,23
+ ,61
+ ,11
+ ,11
+ ,11
+ ,21
+ ,52
+ ,14
+ ,11
+ ,12
+ ,26
+ ,46
+ ,14
+ ,9
+ ,6
+ ,24
+ ,55
+ ,6
+ ,17
+ ,8
+ ,25
+ ,52
+ ,10
+ ,21
+ ,12
+ ,26
+ ,76
+ ,9
+ ,16
+ ,9
+ ,23
+ ,49
+ ,11
+ ,21
+ ,7
+ ,29
+ ,30
+ ,14
+ ,14
+ ,8
+ ,24
+ ,75
+ ,8
+ ,24
+ ,20
+ ,20
+ ,51
+ ,11
+ ,7
+ ,8
+ ,23
+ ,50
+ ,10
+ ,9
+ ,6
+ ,29
+ ,38
+ ,16
+ ,18
+ ,16
+ ,24
+ ,47
+ ,8
+ ,14
+ ,6
+ ,22
+ ,52
+ ,11
+ ,13
+ ,6
+ ,22
+ ,66
+ ,11
+ ,13
+ ,6
+ ,22
+ ,66
+ ,7
+ ,18
+ ,11
+ ,17
+ ,33
+ ,13
+ ,14
+ ,12
+ ,24
+ ,48
+ ,10
+ ,12
+ ,8
+ ,21
+ ,57
+ ,9
+ ,12
+ ,8
+ ,24
+ ,64
+ ,9
+ ,9
+ ,7
+ ,23
+ ,58
+ ,15
+ ,11
+ ,9
+ ,21
+ ,59
+ ,13
+ ,8
+ ,9
+ ,24
+ ,42
+ ,16
+ ,5
+ ,4
+ ,24
+ ,39
+ ,11
+ ,9
+ ,6
+ ,19
+ ,59
+ ,6
+ ,11
+ ,8
+ ,26
+ ,37
+ ,14
+ ,11
+ ,8
+ ,24
+ ,49
+ ,4
+ ,15
+ ,4
+ ,28
+ ,80
+ ,12
+ ,16
+ ,14
+ ,22
+ ,62
+ ,10
+ ,12
+ ,8
+ ,23
+ ,44
+ ,14
+ ,14
+ ,10
+ ,24
+ ,53
+ ,9
+ ,13
+ ,6
+ ,23
+ ,58
+ ,10
+ ,10
+ ,8
+ ,23
+ ,69
+ ,14
+ ,18
+ ,10
+ ,30
+ ,63
+ ,14
+ ,17
+ ,11
+ ,20
+ ,36
+ ,10
+ ,12
+ ,8
+ ,23
+ ,38
+ ,9
+ ,13
+ ,8
+ ,21
+ ,46
+ ,14
+ ,13
+ ,10
+ ,27
+ ,56
+ ,8
+ ,11
+ ,8
+ ,12
+ ,37
+ ,9
+ ,13
+ ,10
+ ,15
+ ,51
+ ,8
+ ,12
+ ,7
+ ,22
+ ,44
+ ,10
+ ,12
+ ,8
+ ,27
+ ,58
+ ,9
+ ,12
+ ,8
+ ,21
+ ,37
+ ,9
+ ,12
+ ,7
+ ,21
+ ,65
+ ,9
+ ,13
+ ,6
+ ,21
+ ,48
+ ,9
+ ,17
+ ,9
+ ,21
+ ,53
+ ,11
+ ,18
+ ,5
+ ,18
+ ,51
+ ,15
+ ,7
+ ,5
+ ,24
+ ,39
+ ,8
+ ,17
+ ,7
+ ,24
+ ,64
+ ,12
+ ,14
+ ,7
+ ,28
+ ,47
+ ,8
+ ,12
+ ,7
+ ,25
+ ,47
+ ,14
+ ,9
+ ,9
+ ,14
+ ,64
+ ,11
+ ,9
+ ,5
+ ,30
+ ,59
+ ,10
+ ,13
+ ,8
+ ,19
+ ,54
+ ,12
+ ,10
+ ,8
+ ,29
+ ,55
+ ,9
+ ,12
+ ,9
+ ,25
+ ,72
+ ,13
+ ,10
+ ,6
+ ,25
+ ,58
+ ,14
+ ,11
+ ,8
+ ,25
+ ,59
+ ,15
+ ,13
+ ,8
+ ,16
+ ,36
+ ,8
+ ,6
+ ,6
+ ,25
+ ,62
+ ,7
+ ,7
+ ,4
+ ,28
+ ,63
+ ,10
+ ,13
+ ,6
+ ,24
+ ,50
+ ,10
+ ,11
+ ,5
+ ,24
+ ,70
+ ,11
+ ,9
+ ,6
+ ,22
+ ,59
+ ,8
+ ,9
+ ,11
+ ,20
+ ,73
+ ,9
+ ,11
+ ,10
+ ,27
+ ,62
+ ,10
+ ,15
+ ,10
+ ,21
+ ,41
+ ,11
+ ,11
+ ,8
+ ,26
+ ,56
+ ,10
+ ,14
+ ,9
+ ,26
+ ,52
+ ,16
+ ,14
+ ,9
+ ,25
+ ,54
+ ,11
+ ,8
+ ,4
+ ,13
+ ,73
+ ,16
+ ,12
+ ,7
+ ,22
+ ,40
+ ,6
+ ,8
+ ,11
+ ,23
+ ,41
+ ,11
+ ,11
+ ,8
+ ,25
+ ,54
+ ,12
+ ,10
+ ,8
+ ,15
+ ,42
+ ,12
+ ,11
+ ,8
+ ,25
+ ,70
+ ,14
+ ,17
+ ,7
+ ,21
+ ,51
+ ,9
+ ,16
+ ,5
+ ,23
+ ,60
+ ,11
+ ,13
+ ,7
+ ,25
+ ,49
+ ,8
+ ,15
+ ,9
+ ,24
+ ,52)
+ ,dim=c(5
+ ,86)
+ ,dimnames=list(c('Doubts'
+ ,'PerantalExpectations'
+ ,'ParentalCriticism'
+ ,'Organization'
+ ,'Intrinsic')
+ ,1:86))
> y <- array(NA,dim=c(5,86),dimnames=list(c('Doubts','PerantalExpectations','ParentalCriticism','Organization','Intrinsic'),1:86))
> 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 = '5'
> #'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
Intrinsic Doubts PerantalExpectations ParentalCriticism Organization
1 68 9 15 6 25
2 48 14 10 8 23
3 44 8 10 7 17
4 67 8 12 9 19
5 46 14 9 8 29
6 54 15 18 11 23
7 61 9 14 9 23
8 52 11 11 11 21
9 46 14 11 12 26
10 55 14 9 6 24
11 52 6 17 8 25
12 76 10 21 12 26
13 49 9 16 9 23
14 30 11 21 7 29
15 75 14 14 8 24
16 51 8 24 20 20
17 50 11 7 8 23
18 38 10 9 6 29
19 47 16 18 16 24
20 52 8 14 6 22
21 66 11 13 6 22
22 66 11 13 6 22
23 33 7 18 11 17
24 48 13 14 12 24
25 57 10 12 8 21
26 64 9 12 8 24
27 58 9 9 7 23
28 59 15 11 9 21
29 42 13 8 9 24
30 39 16 5 4 24
31 59 11 9 6 19
32 37 6 11 8 26
33 49 14 11 8 24
34 80 4 15 4 28
35 62 12 16 14 22
36 44 10 12 8 23
37 53 14 14 10 24
38 58 9 13 6 23
39 69 10 10 8 23
40 63 14 18 10 30
41 36 14 17 11 20
42 38 10 12 8 23
43 46 9 13 8 21
44 56 14 13 10 27
45 37 8 11 8 12
46 51 9 13 10 15
47 44 8 12 7 22
48 58 10 12 8 27
49 37 9 12 8 21
50 65 9 12 7 21
51 48 9 13 6 21
52 53 9 17 9 21
53 51 11 18 5 18
54 39 15 7 5 24
55 64 8 17 7 24
56 47 12 14 7 28
57 47 8 12 7 25
58 64 14 9 9 14
59 59 11 9 5 30
60 54 10 13 8 19
61 55 12 10 8 29
62 72 9 12 9 25
63 58 13 10 6 25
64 59 14 11 8 25
65 36 15 13 8 16
66 62 8 6 6 25
67 63 7 7 4 28
68 50 10 13 6 24
69 70 10 11 5 24
70 59 11 9 6 22
71 73 8 9 11 20
72 62 9 11 10 27
73 41 10 15 10 21
74 56 11 11 8 26
75 52 10 14 9 26
76 54 16 14 9 25
77 73 11 8 4 13
78 40 16 12 7 22
79 41 6 8 11 23
80 54 11 11 8 25
81 42 12 10 8 15
82 70 12 11 8 25
83 51 14 17 7 21
84 60 9 16 5 23
85 49 11 13 7 25
86 52 8 15 9 24
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Doubts PerantalExpectations
54.23866 -0.58774 0.05248
ParentalCriticism Organization
-0.42355 0.37189
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-26.6955 -7.8951 0.2623 6.5570 22.7181
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 54.23866 9.99805 5.425 5.87e-07 ***
Doubts -0.58774 0.45220 -1.300 0.197
PerantalExpectations 0.05248 0.39545 0.133 0.895
ParentalCriticism -0.42355 0.54695 -0.774 0.441
Organization 0.37189 0.32555 1.142 0.257
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 11.01 on 81 degrees of freedom
Multiple R-squared: 0.04836, Adjusted R-squared: 0.00136
F-statistic: 1.029 on 4 and 81 DF, p-value: 0.3975
> 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.2862380 0.5724760 0.71376201
[2,] 0.1681260 0.3362520 0.83187401
[3,] 0.1691051 0.3382102 0.83089488
[4,] 0.2763966 0.5527932 0.72360339
[5,] 0.3254499 0.6508998 0.67455012
[6,] 0.3414471 0.6828942 0.65855290
[7,] 0.7792697 0.4414606 0.22073028
[8,] 0.9178296 0.1643409 0.08217044
[9,] 0.9072519 0.1854963 0.09274813
[10,] 0.8670270 0.2659460 0.13297300
[11,] 0.8696752 0.2606495 0.13032475
[12,] 0.8353864 0.3292273 0.16461365
[13,] 0.7890538 0.4218923 0.21094616
[14,] 0.7685239 0.4629521 0.23147607
[15,] 0.7427557 0.5144885 0.25724425
[16,] 0.8984058 0.2031885 0.10159423
[17,] 0.8645298 0.2709404 0.13547021
[18,] 0.8237446 0.3525109 0.17625543
[19,] 0.8199980 0.3600040 0.18000202
[20,] 0.7773124 0.4453751 0.22268755
[21,] 0.7394904 0.5210192 0.26050960
[22,] 0.7229309 0.5541382 0.27706912
[23,] 0.7642123 0.4715753 0.23578767
[24,] 0.7165251 0.5669499 0.28347494
[25,] 0.7794562 0.4410875 0.22054375
[26,] 0.7322272 0.5355456 0.26777279
[27,] 0.8622760 0.2754480 0.13772400
[28,] 0.8717646 0.2564707 0.12823536
[29,] 0.8640010 0.2719980 0.13599900
[30,] 0.8258384 0.3483231 0.17416156
[31,] 0.7819124 0.4361753 0.21808764
[32,] 0.8231851 0.3536298 0.17681488
[33,] 0.8245668 0.3508663 0.17543316
[34,] 0.8482423 0.3035155 0.15175773
[35,] 0.8838882 0.2322237 0.11611184
[36,] 0.8665132 0.2669736 0.13348682
[37,] 0.8343257 0.3313486 0.16567431
[38,] 0.8676229 0.2647542 0.13237708
[39,] 0.8298678 0.3402644 0.17013221
[40,] 0.8381855 0.3236290 0.16181448
[41,] 0.7976693 0.4046615 0.20233074
[42,] 0.8685645 0.2628709 0.13143546
[43,] 0.8600285 0.2799430 0.13997150
[44,] 0.8453678 0.3092643 0.15463215
[45,] 0.8014229 0.3971542 0.19857710
[46,] 0.7548514 0.4902971 0.24514856
[47,] 0.8165436 0.3669129 0.18345643
[48,] 0.8136906 0.3726188 0.18630938
[49,] 0.7845959 0.4308083 0.21540414
[50,] 0.7802580 0.4394839 0.21974197
[51,] 0.8150815 0.3698371 0.18491854
[52,] 0.7743709 0.4512582 0.22562912
[53,] 0.7155439 0.5689122 0.28445608
[54,] 0.6546357 0.6907286 0.34536431
[55,] 0.7563405 0.4873190 0.24365949
[56,] 0.6943457 0.6113085 0.30565427
[57,] 0.6434713 0.7130574 0.35652870
[58,] 0.6673351 0.6653299 0.33266493
[59,] 0.6026833 0.7946334 0.39731671
[60,] 0.5678379 0.8643242 0.43216211
[61,] 0.5399905 0.9200190 0.46000949
[62,] 0.4929851 0.9859703 0.50701486
[63,] 0.4151524 0.8303047 0.58484764
[64,] 0.7635431 0.4729138 0.23645689
[65,] 0.7679140 0.4641719 0.23208595
[66,] 0.6854845 0.6290311 0.31451553
[67,] 0.5758423 0.8483153 0.42415766
[68,] 0.4592462 0.9184924 0.54075382
[69,] 0.4396618 0.8793237 0.56033817
[70,] 0.4624351 0.9248702 0.53756492
[71,] 0.5381714 0.9236572 0.46182858
> postscript(file="/var/www/html/freestat/rcomp/tmp/1hy6i1290270350.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/freestat/rcomp/tmp/2a7n31290270350.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/freestat/rcomp/tmp/3a7n31290270350.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/freestat/rcomp/tmp/4a7n31290270350.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/freestat/rcomp/tmp/5lh4o1290270350.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 = 86
Frequency = 1
1 2 3 4 5 6
11.5078962 -3.7001218 -9.4187843 13.5795793 -7.8790011 3.7384699
7 8 9 10 11 12
6.5748041 0.4986021 -5.1740769 2.1333623 -5.5131892 21.9501936
13 14 15 16 17 18
-5.5301461 -26.6954856 22.7180854 1.2370281 -3.3059303 -19.0770782
19 20 21 22 23 24
-0.9279313 -3.9116961 11.9040129 11.9040129 -19.7321325 -3.1754620
25 26 27 28 29 30
3.5877344 8.8843126 2.9900810 9.0024819 -10.1312593 -13.3283467
31 32 33 34 35 36
6.2295905 -20.5702310 -3.1244893 18.6063973 11.7227267 -10.1560504
37 38 39 40 41 42
1.5651840 2.3566312 14.9488998 9.1239291 -13.6811223 -16.1560504
43 44 45 46 47 48
-8.0524853 3.5019819 -14.1882480 0.0259678 -11.3831965 2.3563799
49 50 51 52 53 54
-17.0000102 10.5764405 -6.8995839 -0.8388364 -2.2943423 -13.5974922
55 56 57 58 59 60
7.6106431 -8.3685228 -9.4988738 16.1229344 1.7152247 1.2790441
61 62 63 64 65 66
-0.1069655 16.9359695 4.1212502 6.5036183 -12.6665556 5.3924275
67 68 69 70 71 72
3.7894320 -5.4275166 14.2538843 5.1139133 20.2122109 6.6682090
73 74 75 76 77 78
-11.7225923 1.3684921 -2.9531286 2.9452315 21.6663215 -10.6812397
79 80 81 82 83 84
-14.0264805 -0.2596155 -7.9004717 16.3281291 -0.7472120 3.7756566
85 86
-5.7881151 -3.4373081
> postscript(file="/var/www/html/freestat/rcomp/tmp/6lh4o1290270350.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 = 86
Frequency = 1
lag(myerror, k = 1) myerror
0 11.5078962 NA
1 -3.7001218 11.5078962
2 -9.4187843 -3.7001218
3 13.5795793 -9.4187843
4 -7.8790011 13.5795793
5 3.7384699 -7.8790011
6 6.5748041 3.7384699
7 0.4986021 6.5748041
8 -5.1740769 0.4986021
9 2.1333623 -5.1740769
10 -5.5131892 2.1333623
11 21.9501936 -5.5131892
12 -5.5301461 21.9501936
13 -26.6954856 -5.5301461
14 22.7180854 -26.6954856
15 1.2370281 22.7180854
16 -3.3059303 1.2370281
17 -19.0770782 -3.3059303
18 -0.9279313 -19.0770782
19 -3.9116961 -0.9279313
20 11.9040129 -3.9116961
21 11.9040129 11.9040129
22 -19.7321325 11.9040129
23 -3.1754620 -19.7321325
24 3.5877344 -3.1754620
25 8.8843126 3.5877344
26 2.9900810 8.8843126
27 9.0024819 2.9900810
28 -10.1312593 9.0024819
29 -13.3283467 -10.1312593
30 6.2295905 -13.3283467
31 -20.5702310 6.2295905
32 -3.1244893 -20.5702310
33 18.6063973 -3.1244893
34 11.7227267 18.6063973
35 -10.1560504 11.7227267
36 1.5651840 -10.1560504
37 2.3566312 1.5651840
38 14.9488998 2.3566312
39 9.1239291 14.9488998
40 -13.6811223 9.1239291
41 -16.1560504 -13.6811223
42 -8.0524853 -16.1560504
43 3.5019819 -8.0524853
44 -14.1882480 3.5019819
45 0.0259678 -14.1882480
46 -11.3831965 0.0259678
47 2.3563799 -11.3831965
48 -17.0000102 2.3563799
49 10.5764405 -17.0000102
50 -6.8995839 10.5764405
51 -0.8388364 -6.8995839
52 -2.2943423 -0.8388364
53 -13.5974922 -2.2943423
54 7.6106431 -13.5974922
55 -8.3685228 7.6106431
56 -9.4988738 -8.3685228
57 16.1229344 -9.4988738
58 1.7152247 16.1229344
59 1.2790441 1.7152247
60 -0.1069655 1.2790441
61 16.9359695 -0.1069655
62 4.1212502 16.9359695
63 6.5036183 4.1212502
64 -12.6665556 6.5036183
65 5.3924275 -12.6665556
66 3.7894320 5.3924275
67 -5.4275166 3.7894320
68 14.2538843 -5.4275166
69 5.1139133 14.2538843
70 20.2122109 5.1139133
71 6.6682090 20.2122109
72 -11.7225923 6.6682090
73 1.3684921 -11.7225923
74 -2.9531286 1.3684921
75 2.9452315 -2.9531286
76 21.6663215 2.9452315
77 -10.6812397 21.6663215
78 -14.0264805 -10.6812397
79 -0.2596155 -14.0264805
80 -7.9004717 -0.2596155
81 16.3281291 -7.9004717
82 -0.7472120 16.3281291
83 3.7756566 -0.7472120
84 -5.7881151 3.7756566
85 -3.4373081 -5.7881151
86 NA -3.4373081
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -3.7001218 11.5078962
[2,] -9.4187843 -3.7001218
[3,] 13.5795793 -9.4187843
[4,] -7.8790011 13.5795793
[5,] 3.7384699 -7.8790011
[6,] 6.5748041 3.7384699
[7,] 0.4986021 6.5748041
[8,] -5.1740769 0.4986021
[9,] 2.1333623 -5.1740769
[10,] -5.5131892 2.1333623
[11,] 21.9501936 -5.5131892
[12,] -5.5301461 21.9501936
[13,] -26.6954856 -5.5301461
[14,] 22.7180854 -26.6954856
[15,] 1.2370281 22.7180854
[16,] -3.3059303 1.2370281
[17,] -19.0770782 -3.3059303
[18,] -0.9279313 -19.0770782
[19,] -3.9116961 -0.9279313
[20,] 11.9040129 -3.9116961
[21,] 11.9040129 11.9040129
[22,] -19.7321325 11.9040129
[23,] -3.1754620 -19.7321325
[24,] 3.5877344 -3.1754620
[25,] 8.8843126 3.5877344
[26,] 2.9900810 8.8843126
[27,] 9.0024819 2.9900810
[28,] -10.1312593 9.0024819
[29,] -13.3283467 -10.1312593
[30,] 6.2295905 -13.3283467
[31,] -20.5702310 6.2295905
[32,] -3.1244893 -20.5702310
[33,] 18.6063973 -3.1244893
[34,] 11.7227267 18.6063973
[35,] -10.1560504 11.7227267
[36,] 1.5651840 -10.1560504
[37,] 2.3566312 1.5651840
[38,] 14.9488998 2.3566312
[39,] 9.1239291 14.9488998
[40,] -13.6811223 9.1239291
[41,] -16.1560504 -13.6811223
[42,] -8.0524853 -16.1560504
[43,] 3.5019819 -8.0524853
[44,] -14.1882480 3.5019819
[45,] 0.0259678 -14.1882480
[46,] -11.3831965 0.0259678
[47,] 2.3563799 -11.3831965
[48,] -17.0000102 2.3563799
[49,] 10.5764405 -17.0000102
[50,] -6.8995839 10.5764405
[51,] -0.8388364 -6.8995839
[52,] -2.2943423 -0.8388364
[53,] -13.5974922 -2.2943423
[54,] 7.6106431 -13.5974922
[55,] -8.3685228 7.6106431
[56,] -9.4988738 -8.3685228
[57,] 16.1229344 -9.4988738
[58,] 1.7152247 16.1229344
[59,] 1.2790441 1.7152247
[60,] -0.1069655 1.2790441
[61,] 16.9359695 -0.1069655
[62,] 4.1212502 16.9359695
[63,] 6.5036183 4.1212502
[64,] -12.6665556 6.5036183
[65,] 5.3924275 -12.6665556
[66,] 3.7894320 5.3924275
[67,] -5.4275166 3.7894320
[68,] 14.2538843 -5.4275166
[69,] 5.1139133 14.2538843
[70,] 20.2122109 5.1139133
[71,] 6.6682090 20.2122109
[72,] -11.7225923 6.6682090
[73,] 1.3684921 -11.7225923
[74,] -2.9531286 1.3684921
[75,] 2.9452315 -2.9531286
[76,] 21.6663215 2.9452315
[77,] -10.6812397 21.6663215
[78,] -14.0264805 -10.6812397
[79,] -0.2596155 -14.0264805
[80,] -7.9004717 -0.2596155
[81,] 16.3281291 -7.9004717
[82,] -0.7472120 16.3281291
[83,] 3.7756566 -0.7472120
[84,] -5.7881151 3.7756566
[85,] -3.4373081 -5.7881151
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -3.7001218 11.5078962
2 -9.4187843 -3.7001218
3 13.5795793 -9.4187843
4 -7.8790011 13.5795793
5 3.7384699 -7.8790011
6 6.5748041 3.7384699
7 0.4986021 6.5748041
8 -5.1740769 0.4986021
9 2.1333623 -5.1740769
10 -5.5131892 2.1333623
11 21.9501936 -5.5131892
12 -5.5301461 21.9501936
13 -26.6954856 -5.5301461
14 22.7180854 -26.6954856
15 1.2370281 22.7180854
16 -3.3059303 1.2370281
17 -19.0770782 -3.3059303
18 -0.9279313 -19.0770782
19 -3.9116961 -0.9279313
20 11.9040129 -3.9116961
21 11.9040129 11.9040129
22 -19.7321325 11.9040129
23 -3.1754620 -19.7321325
24 3.5877344 -3.1754620
25 8.8843126 3.5877344
26 2.9900810 8.8843126
27 9.0024819 2.9900810
28 -10.1312593 9.0024819
29 -13.3283467 -10.1312593
30 6.2295905 -13.3283467
31 -20.5702310 6.2295905
32 -3.1244893 -20.5702310
33 18.6063973 -3.1244893
34 11.7227267 18.6063973
35 -10.1560504 11.7227267
36 1.5651840 -10.1560504
37 2.3566312 1.5651840
38 14.9488998 2.3566312
39 9.1239291 14.9488998
40 -13.6811223 9.1239291
41 -16.1560504 -13.6811223
42 -8.0524853 -16.1560504
43 3.5019819 -8.0524853
44 -14.1882480 3.5019819
45 0.0259678 -14.1882480
46 -11.3831965 0.0259678
47 2.3563799 -11.3831965
48 -17.0000102 2.3563799
49 10.5764405 -17.0000102
50 -6.8995839 10.5764405
51 -0.8388364 -6.8995839
52 -2.2943423 -0.8388364
53 -13.5974922 -2.2943423
54 7.6106431 -13.5974922
55 -8.3685228 7.6106431
56 -9.4988738 -8.3685228
57 16.1229344 -9.4988738
58 1.7152247 16.1229344
59 1.2790441 1.7152247
60 -0.1069655 1.2790441
61 16.9359695 -0.1069655
62 4.1212502 16.9359695
63 6.5036183 4.1212502
64 -12.6665556 6.5036183
65 5.3924275 -12.6665556
66 3.7894320 5.3924275
67 -5.4275166 3.7894320
68 14.2538843 -5.4275166
69 5.1139133 14.2538843
70 20.2122109 5.1139133
71 6.6682090 20.2122109
72 -11.7225923 6.6682090
73 1.3684921 -11.7225923
74 -2.9531286 1.3684921
75 2.9452315 -2.9531286
76 21.6663215 2.9452315
77 -10.6812397 21.6663215
78 -14.0264805 -10.6812397
79 -0.2596155 -14.0264805
80 -7.9004717 -0.2596155
81 16.3281291 -7.9004717
82 -0.7472120 16.3281291
83 3.7756566 -0.7472120
84 -5.7881151 3.7756566
85 -3.4373081 -5.7881151
> 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/freestat/rcomp/tmp/7v84r1290270350.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/freestat/rcomp/tmp/8v84r1290270350.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/freestat/rcomp/tmp/96hlu1290270350.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/freestat/rcomp/tmp/106hlu1290270350.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11r0j01290270350.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/freestat/rcomp/tmp/12d0i61290270350.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/freestat/rcomp/tmp/139syw1290270350.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/freestat/rcomp/tmp/14uae21290270350.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/freestat/rcomp/tmp/15ftu81290270350.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/freestat/rcomp/tmp/16t2az1290270350.tab")
+ }
>
> try(system("convert tmp/1hy6i1290270350.ps tmp/1hy6i1290270350.png",intern=TRUE))
character(0)
> try(system("convert tmp/2a7n31290270350.ps tmp/2a7n31290270350.png",intern=TRUE))
character(0)
> try(system("convert tmp/3a7n31290270350.ps tmp/3a7n31290270350.png",intern=TRUE))
character(0)
> try(system("convert tmp/4a7n31290270350.ps tmp/4a7n31290270350.png",intern=TRUE))
character(0)
> try(system("convert tmp/5lh4o1290270350.ps tmp/5lh4o1290270350.png",intern=TRUE))
character(0)
> try(system("convert tmp/6lh4o1290270350.ps tmp/6lh4o1290270350.png",intern=TRUE))
character(0)
> try(system("convert tmp/7v84r1290270350.ps tmp/7v84r1290270350.png",intern=TRUE))
character(0)
> try(system("convert tmp/8v84r1290270350.ps tmp/8v84r1290270350.png",intern=TRUE))
character(0)
> try(system("convert tmp/96hlu1290270350.ps tmp/96hlu1290270350.png",intern=TRUE))
character(0)
> try(system("convert tmp/106hlu1290270350.ps tmp/106hlu1290270350.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.256 2.550 5.326