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(55,1,20,0,80,0,52,0,75,1,30,1,90,1,68,1,24,0,60,0,65,1,60,0,80,1,65,1,90,0,65,1,76,1,70,1,38,0,60,1,10,0,5,0,93,1,70,0,61,1,72,1,40,0,75,1,100,1,29,0,70,1,25,0,70,1,82,0,40,0,50,1,70,1,91,1,10,0,25,0,56,0,30,0,74,0,60,0,80,0,80,1,60,1,64,1,40,1,80,1,71,1,65,1,90,0,68,1,76,1,25,1,79,1,40,0,61,1,27,1,70,0,40,0,13,0,15,0,38,1,47,0,65,1,62,1,50,0,80,1,87,1,40,1,80,1,20,0,60,1,48,1,70,1,91,1,10,0,50,0,70,0,45,1,20,1,10,0,90,1,80,1,74,0,71,0,40,0,29,1,60,1,31,0,67,0,82,0,40,1,30,1,70,0,63,0,35,0,35,1,70,0,60,1,80,1,70,1,71,1),dim=c(2,105),dimnames=list(c('Q1','gender'),1:105))
> y <- array(NA,dim=c(2,105),dimnames=list(c('Q1','gender'),1:105))
> 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 = '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.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
Q1 gender t
1 55 1 1
2 20 0 2
3 80 0 3
4 52 0 4
5 75 1 5
6 30 1 6
7 90 1 7
8 68 1 8
9 24 0 9
10 60 0 10
11 65 1 11
12 60 0 12
13 80 1 13
14 65 1 14
15 90 0 15
16 65 1 16
17 76 1 17
18 70 1 18
19 38 0 19
20 60 1 20
21 10 0 21
22 5 0 22
23 93 1 23
24 70 0 24
25 61 1 25
26 72 1 26
27 40 0 27
28 75 1 28
29 100 1 29
30 29 0 30
31 70 1 31
32 25 0 32
33 70 1 33
34 82 0 34
35 40 0 35
36 50 1 36
37 70 1 37
38 91 1 38
39 10 0 39
40 25 0 40
41 56 0 41
42 30 0 42
43 74 0 43
44 60 0 44
45 80 0 45
46 80 1 46
47 60 1 47
48 64 1 48
49 40 1 49
50 80 1 50
51 71 1 51
52 65 1 52
53 90 0 53
54 68 1 54
55 76 1 55
56 25 1 56
57 79 1 57
58 40 0 58
59 61 1 59
60 27 1 60
61 70 0 61
62 40 0 62
63 13 0 63
64 15 0 64
65 38 1 65
66 47 0 66
67 65 1 67
68 62 1 68
69 50 0 69
70 80 1 70
71 87 1 71
72 40 1 72
73 80 1 73
74 20 0 74
75 60 1 75
76 48 1 76
77 70 1 77
78 91 1 78
79 10 0 79
80 50 0 80
81 70 0 81
82 45 1 82
83 20 1 83
84 10 0 84
85 90 1 85
86 80 1 86
87 74 0 87
88 71 0 88
89 40 0 89
90 29 1 90
91 60 1 91
92 31 0 92
93 67 0 93
94 82 0 94
95 40 1 95
96 30 1 96
97 70 0 97
98 63 0 98
99 35 0 99
100 35 1 100
101 70 0 101
102 60 1 102
103 80 1 103
104 70 1 104
105 71 1 105
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) gender t
50.03134 16.49221 -0.04678
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-44.002 -14.968 1.851 16.751 42.448
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 50.03134 4.91902 10.171 < 2e-16 ***
gender 16.49221 4.30341 3.832 0.000220 ***
t -0.04678 0.07044 -0.664 0.508147
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 21.88 on 102 degrees of freedom
Multiple R-squared: 0.1289, Adjusted R-squared: 0.1118
F-statistic: 7.544 on 2 and 102 DF, p-value: 0.0008803
> 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.8787958 0.2424084 0.12120420
[2,] 0.8757827 0.2484345 0.12421727
[3,] 0.7967264 0.4065472 0.20327360
[4,] 0.8246228 0.3507545 0.17537723
[5,] 0.7675834 0.4648332 0.23241659
[6,] 0.6780790 0.6438419 0.32192097
[7,] 0.5921123 0.8157754 0.40788768
[8,] 0.5099203 0.9801594 0.49007969
[9,] 0.4252777 0.8505554 0.57472229
[10,] 0.4908356 0.9816712 0.50916440
[11,] 0.4303284 0.8606569 0.56967156
[12,] 0.3508507 0.7017014 0.64914929
[13,] 0.2852788 0.5705576 0.71472122
[14,] 0.3124818 0.6249637 0.68751817
[15,] 0.2657219 0.5314438 0.73427810
[16,] 0.4620452 0.9240903 0.53795483
[17,] 0.6090575 0.7818849 0.39094247
[18,] 0.6558849 0.6882303 0.34411515
[19,] 0.6676185 0.6647631 0.33238153
[20,] 0.6084027 0.7831947 0.39159733
[21,] 0.5447922 0.9104156 0.45520781
[22,] 0.4838754 0.9677509 0.51612456
[23,] 0.4258784 0.8517567 0.57412165
[24,] 0.4901355 0.9802709 0.50986454
[25,] 0.4720297 0.9440594 0.52797032
[26,] 0.4109057 0.8218113 0.58909433
[27,] 0.4048238 0.8096476 0.59517621
[28,] 0.3467612 0.6935224 0.65323881
[29,] 0.4451415 0.8902830 0.55485852
[30,] 0.3927787 0.7855575 0.60722127
[31,] 0.3767372 0.7534744 0.62326278
[32,] 0.3227154 0.6454308 0.67728460
[33,] 0.3335558 0.6671116 0.66644420
[34,] 0.4347312 0.8694625 0.56526877
[35,] 0.4285809 0.8571618 0.57141909
[36,] 0.3890185 0.7780370 0.61098148
[37,] 0.3667136 0.7334272 0.63328640
[38,] 0.4037374 0.8074748 0.59626262
[39,] 0.3672814 0.7345628 0.63271858
[40,] 0.4277053 0.8554106 0.57229469
[41,] 0.3943091 0.7886183 0.60569087
[42,] 0.3530198 0.7060396 0.64698021
[43,] 0.3073172 0.6146344 0.69268280
[44,] 0.3301938 0.6603876 0.66980622
[45,] 0.3061098 0.6122196 0.69389020
[46,] 0.2656495 0.5312991 0.73435046
[47,] 0.2249532 0.4499063 0.77504683
[48,] 0.3700888 0.7401777 0.62991116
[49,] 0.3299585 0.6599170 0.67004152
[50,] 0.3109072 0.6218144 0.68909281
[51,] 0.4142638 0.8285276 0.58573621
[52,] 0.4071800 0.8143600 0.59282001
[53,] 0.3565582 0.7131163 0.64344185
[54,] 0.3131526 0.6263051 0.68684743
[55,] 0.3776818 0.7553636 0.62231822
[56,] 0.4035806 0.8071613 0.59641936
[57,] 0.3513172 0.7026343 0.64868283
[58,] 0.3994032 0.7988064 0.60059682
[59,] 0.4441216 0.8882432 0.55587838
[60,] 0.4425110 0.8850220 0.55748901
[61,] 0.3868011 0.7736023 0.61319886
[62,] 0.3339162 0.6678324 0.66608380
[63,] 0.2820327 0.5640655 0.71796727
[64,] 0.2362749 0.4725498 0.76372511
[65,] 0.2301784 0.4603569 0.76982156
[66,] 0.2681341 0.5362683 0.73186586
[67,] 0.2484054 0.4968107 0.75159464
[68,] 0.2588922 0.5177845 0.74110776
[69,] 0.2709133 0.5418267 0.72908666
[70,] 0.2248214 0.4496429 0.77517857
[71,] 0.1864267 0.3728535 0.81357326
[72,] 0.1650580 0.3301160 0.83494199
[73,] 0.2786889 0.5573777 0.72131114
[74,] 0.3686412 0.7372824 0.63135881
[75,] 0.3102255 0.6204511 0.68977445
[76,] 0.3198876 0.6397752 0.68011240
[77,] 0.2686973 0.5373945 0.73130275
[78,] 0.3375066 0.6750132 0.66249341
[79,] 0.5566015 0.8867969 0.44339847
[80,] 0.6694412 0.6611175 0.33055877
[81,] 0.7799729 0.4400541 0.22002707
[82,] 0.8183887 0.3632227 0.18161134
[83,] 0.8576006 0.2847987 0.14239937
[84,] 0.8046089 0.3907823 0.19539114
[85,] 0.7658254 0.4683493 0.23417464
[86,] 0.7904154 0.4191692 0.20958459
[87,] 0.7850157 0.4299686 0.21498432
[88,] 0.7467948 0.5064104 0.25320520
[89,] 0.9042903 0.1914193 0.09570966
[90,] 0.8655015 0.2689969 0.13449845
[91,] 0.7893795 0.4212409 0.21062045
[92,] 0.8456794 0.3086412 0.15432060
[93,] 0.8903105 0.2193790 0.10968952
[94,] 0.8575190 0.2849620 0.14248100
> postscript(file="/var/www/html/freestat/rcomp/tmp/11igb1290508399.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/21igb1290508399.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/3uafe1290508399.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/4uafe1290508399.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/5uafe1290508399.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 = 105
Frequency = 1
1 2 3 4 5 6
-11.47676572 -29.93777691 30.10900309 2.15578309 8.71035428 -36.24286572
7 8 9 10 11 12
23.80391428 1.85069428 -25.61031692 10.43646308 -1.00896573 10.53002308
13 14 15 16 17 18
14.08459427 -0.86862573 40.67036308 -0.77506573 10.27171427 4.31849427
19 20 21 22 23 24
-11.14251693 -5.58794573 -39.04895693 -44.00217693 27.55239426 21.09138307
25 26 27 28 29 30
-4.35404574 6.69273426 -8.76827693 9.78629426 34.83307426 -19.62793694
31 32 33 34 35 36
4.92663426 -23.53437694 5.02019426 33.55918306 -8.39403694 -14.83946575
37 38 39 40 41 42
5.20731425 26.25409425 -38.20691694 -23.16013695 7.88664305 -18.06657695
43 44 45 46 47 48
25.98020305 12.02698305 32.07376305 15.62833424 -4.32488576 -0.27810576
49 50 51 52 53 54
-24.23132576 15.81545424 6.86223424 0.90901424 42.44800304 4.00257424
55 56 57 58 59 60
12.04935424 -38.90386577 15.14291423 -7.31809696 -2.76352577 -36.71674577
61 62 63 64 65 66
22.82224304 -7.13097696 -34.08419697 -32.03741697 -25.48284577 0.05614303
67 68 69 70 71 72
1.61071423 -1.34250578 3.19648303 16.75105422 23.79783422 -23.15538578
73 74 75 76 77 78
16.89139422 -26.56961698 -3.01504578 -14.96826578 7.07851422 28.12529422
79 80 81 82 83 84
-36.33571698 3.71106302 23.75784302 -17.68758579 -42.64080579 -36.10181698
85 86 87 88 89 90
27.45275421 17.49953421 28.03852301 25.08530301 -5.86791699 -33.31334580
91 92 93 94 95 96
-2.26656580 -14.72757699 21.31920301 36.36598301 -22.07944580 -32.03266580
97 98 99 100 101 102
24.50632300 17.55310300 -10.40011700 -26.84554580 24.69344300 -1.75198581
103 104 105
18.29479419 8.34157419 9.38835419
> postscript(file="/var/www/html/freestat/rcomp/tmp/6mjfy1290508399.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 = 105
Frequency = 1
lag(myerror, k = 1) myerror
0 -11.47676572 NA
1 -29.93777691 -11.47676572
2 30.10900309 -29.93777691
3 2.15578309 30.10900309
4 8.71035428 2.15578309
5 -36.24286572 8.71035428
6 23.80391428 -36.24286572
7 1.85069428 23.80391428
8 -25.61031692 1.85069428
9 10.43646308 -25.61031692
10 -1.00896573 10.43646308
11 10.53002308 -1.00896573
12 14.08459427 10.53002308
13 -0.86862573 14.08459427
14 40.67036308 -0.86862573
15 -0.77506573 40.67036308
16 10.27171427 -0.77506573
17 4.31849427 10.27171427
18 -11.14251693 4.31849427
19 -5.58794573 -11.14251693
20 -39.04895693 -5.58794573
21 -44.00217693 -39.04895693
22 27.55239426 -44.00217693
23 21.09138307 27.55239426
24 -4.35404574 21.09138307
25 6.69273426 -4.35404574
26 -8.76827693 6.69273426
27 9.78629426 -8.76827693
28 34.83307426 9.78629426
29 -19.62793694 34.83307426
30 4.92663426 -19.62793694
31 -23.53437694 4.92663426
32 5.02019426 -23.53437694
33 33.55918306 5.02019426
34 -8.39403694 33.55918306
35 -14.83946575 -8.39403694
36 5.20731425 -14.83946575
37 26.25409425 5.20731425
38 -38.20691694 26.25409425
39 -23.16013695 -38.20691694
40 7.88664305 -23.16013695
41 -18.06657695 7.88664305
42 25.98020305 -18.06657695
43 12.02698305 25.98020305
44 32.07376305 12.02698305
45 15.62833424 32.07376305
46 -4.32488576 15.62833424
47 -0.27810576 -4.32488576
48 -24.23132576 -0.27810576
49 15.81545424 -24.23132576
50 6.86223424 15.81545424
51 0.90901424 6.86223424
52 42.44800304 0.90901424
53 4.00257424 42.44800304
54 12.04935424 4.00257424
55 -38.90386577 12.04935424
56 15.14291423 -38.90386577
57 -7.31809696 15.14291423
58 -2.76352577 -7.31809696
59 -36.71674577 -2.76352577
60 22.82224304 -36.71674577
61 -7.13097696 22.82224304
62 -34.08419697 -7.13097696
63 -32.03741697 -34.08419697
64 -25.48284577 -32.03741697
65 0.05614303 -25.48284577
66 1.61071423 0.05614303
67 -1.34250578 1.61071423
68 3.19648303 -1.34250578
69 16.75105422 3.19648303
70 23.79783422 16.75105422
71 -23.15538578 23.79783422
72 16.89139422 -23.15538578
73 -26.56961698 16.89139422
74 -3.01504578 -26.56961698
75 -14.96826578 -3.01504578
76 7.07851422 -14.96826578
77 28.12529422 7.07851422
78 -36.33571698 28.12529422
79 3.71106302 -36.33571698
80 23.75784302 3.71106302
81 -17.68758579 23.75784302
82 -42.64080579 -17.68758579
83 -36.10181698 -42.64080579
84 27.45275421 -36.10181698
85 17.49953421 27.45275421
86 28.03852301 17.49953421
87 25.08530301 28.03852301
88 -5.86791699 25.08530301
89 -33.31334580 -5.86791699
90 -2.26656580 -33.31334580
91 -14.72757699 -2.26656580
92 21.31920301 -14.72757699
93 36.36598301 21.31920301
94 -22.07944580 36.36598301
95 -32.03266580 -22.07944580
96 24.50632300 -32.03266580
97 17.55310300 24.50632300
98 -10.40011700 17.55310300
99 -26.84554580 -10.40011700
100 24.69344300 -26.84554580
101 -1.75198581 24.69344300
102 18.29479419 -1.75198581
103 8.34157419 18.29479419
104 9.38835419 8.34157419
105 NA 9.38835419
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -29.93777691 -11.47676572
[2,] 30.10900309 -29.93777691
[3,] 2.15578309 30.10900309
[4,] 8.71035428 2.15578309
[5,] -36.24286572 8.71035428
[6,] 23.80391428 -36.24286572
[7,] 1.85069428 23.80391428
[8,] -25.61031692 1.85069428
[9,] 10.43646308 -25.61031692
[10,] -1.00896573 10.43646308
[11,] 10.53002308 -1.00896573
[12,] 14.08459427 10.53002308
[13,] -0.86862573 14.08459427
[14,] 40.67036308 -0.86862573
[15,] -0.77506573 40.67036308
[16,] 10.27171427 -0.77506573
[17,] 4.31849427 10.27171427
[18,] -11.14251693 4.31849427
[19,] -5.58794573 -11.14251693
[20,] -39.04895693 -5.58794573
[21,] -44.00217693 -39.04895693
[22,] 27.55239426 -44.00217693
[23,] 21.09138307 27.55239426
[24,] -4.35404574 21.09138307
[25,] 6.69273426 -4.35404574
[26,] -8.76827693 6.69273426
[27,] 9.78629426 -8.76827693
[28,] 34.83307426 9.78629426
[29,] -19.62793694 34.83307426
[30,] 4.92663426 -19.62793694
[31,] -23.53437694 4.92663426
[32,] 5.02019426 -23.53437694
[33,] 33.55918306 5.02019426
[34,] -8.39403694 33.55918306
[35,] -14.83946575 -8.39403694
[36,] 5.20731425 -14.83946575
[37,] 26.25409425 5.20731425
[38,] -38.20691694 26.25409425
[39,] -23.16013695 -38.20691694
[40,] 7.88664305 -23.16013695
[41,] -18.06657695 7.88664305
[42,] 25.98020305 -18.06657695
[43,] 12.02698305 25.98020305
[44,] 32.07376305 12.02698305
[45,] 15.62833424 32.07376305
[46,] -4.32488576 15.62833424
[47,] -0.27810576 -4.32488576
[48,] -24.23132576 -0.27810576
[49,] 15.81545424 -24.23132576
[50,] 6.86223424 15.81545424
[51,] 0.90901424 6.86223424
[52,] 42.44800304 0.90901424
[53,] 4.00257424 42.44800304
[54,] 12.04935424 4.00257424
[55,] -38.90386577 12.04935424
[56,] 15.14291423 -38.90386577
[57,] -7.31809696 15.14291423
[58,] -2.76352577 -7.31809696
[59,] -36.71674577 -2.76352577
[60,] 22.82224304 -36.71674577
[61,] -7.13097696 22.82224304
[62,] -34.08419697 -7.13097696
[63,] -32.03741697 -34.08419697
[64,] -25.48284577 -32.03741697
[65,] 0.05614303 -25.48284577
[66,] 1.61071423 0.05614303
[67,] -1.34250578 1.61071423
[68,] 3.19648303 -1.34250578
[69,] 16.75105422 3.19648303
[70,] 23.79783422 16.75105422
[71,] -23.15538578 23.79783422
[72,] 16.89139422 -23.15538578
[73,] -26.56961698 16.89139422
[74,] -3.01504578 -26.56961698
[75,] -14.96826578 -3.01504578
[76,] 7.07851422 -14.96826578
[77,] 28.12529422 7.07851422
[78,] -36.33571698 28.12529422
[79,] 3.71106302 -36.33571698
[80,] 23.75784302 3.71106302
[81,] -17.68758579 23.75784302
[82,] -42.64080579 -17.68758579
[83,] -36.10181698 -42.64080579
[84,] 27.45275421 -36.10181698
[85,] 17.49953421 27.45275421
[86,] 28.03852301 17.49953421
[87,] 25.08530301 28.03852301
[88,] -5.86791699 25.08530301
[89,] -33.31334580 -5.86791699
[90,] -2.26656580 -33.31334580
[91,] -14.72757699 -2.26656580
[92,] 21.31920301 -14.72757699
[93,] 36.36598301 21.31920301
[94,] -22.07944580 36.36598301
[95,] -32.03266580 -22.07944580
[96,] 24.50632300 -32.03266580
[97,] 17.55310300 24.50632300
[98,] -10.40011700 17.55310300
[99,] -26.84554580 -10.40011700
[100,] 24.69344300 -26.84554580
[101,] -1.75198581 24.69344300
[102,] 18.29479419 -1.75198581
[103,] 8.34157419 18.29479419
[104,] 9.38835419 8.34157419
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -29.93777691 -11.47676572
2 30.10900309 -29.93777691
3 2.15578309 30.10900309
4 8.71035428 2.15578309
5 -36.24286572 8.71035428
6 23.80391428 -36.24286572
7 1.85069428 23.80391428
8 -25.61031692 1.85069428
9 10.43646308 -25.61031692
10 -1.00896573 10.43646308
11 10.53002308 -1.00896573
12 14.08459427 10.53002308
13 -0.86862573 14.08459427
14 40.67036308 -0.86862573
15 -0.77506573 40.67036308
16 10.27171427 -0.77506573
17 4.31849427 10.27171427
18 -11.14251693 4.31849427
19 -5.58794573 -11.14251693
20 -39.04895693 -5.58794573
21 -44.00217693 -39.04895693
22 27.55239426 -44.00217693
23 21.09138307 27.55239426
24 -4.35404574 21.09138307
25 6.69273426 -4.35404574
26 -8.76827693 6.69273426
27 9.78629426 -8.76827693
28 34.83307426 9.78629426
29 -19.62793694 34.83307426
30 4.92663426 -19.62793694
31 -23.53437694 4.92663426
32 5.02019426 -23.53437694
33 33.55918306 5.02019426
34 -8.39403694 33.55918306
35 -14.83946575 -8.39403694
36 5.20731425 -14.83946575
37 26.25409425 5.20731425
38 -38.20691694 26.25409425
39 -23.16013695 -38.20691694
40 7.88664305 -23.16013695
41 -18.06657695 7.88664305
42 25.98020305 -18.06657695
43 12.02698305 25.98020305
44 32.07376305 12.02698305
45 15.62833424 32.07376305
46 -4.32488576 15.62833424
47 -0.27810576 -4.32488576
48 -24.23132576 -0.27810576
49 15.81545424 -24.23132576
50 6.86223424 15.81545424
51 0.90901424 6.86223424
52 42.44800304 0.90901424
53 4.00257424 42.44800304
54 12.04935424 4.00257424
55 -38.90386577 12.04935424
56 15.14291423 -38.90386577
57 -7.31809696 15.14291423
58 -2.76352577 -7.31809696
59 -36.71674577 -2.76352577
60 22.82224304 -36.71674577
61 -7.13097696 22.82224304
62 -34.08419697 -7.13097696
63 -32.03741697 -34.08419697
64 -25.48284577 -32.03741697
65 0.05614303 -25.48284577
66 1.61071423 0.05614303
67 -1.34250578 1.61071423
68 3.19648303 -1.34250578
69 16.75105422 3.19648303
70 23.79783422 16.75105422
71 -23.15538578 23.79783422
72 16.89139422 -23.15538578
73 -26.56961698 16.89139422
74 -3.01504578 -26.56961698
75 -14.96826578 -3.01504578
76 7.07851422 -14.96826578
77 28.12529422 7.07851422
78 -36.33571698 28.12529422
79 3.71106302 -36.33571698
80 23.75784302 3.71106302
81 -17.68758579 23.75784302
82 -42.64080579 -17.68758579
83 -36.10181698 -42.64080579
84 27.45275421 -36.10181698
85 17.49953421 27.45275421
86 28.03852301 17.49953421
87 25.08530301 28.03852301
88 -5.86791699 25.08530301
89 -33.31334580 -5.86791699
90 -2.26656580 -33.31334580
91 -14.72757699 -2.26656580
92 21.31920301 -14.72757699
93 36.36598301 21.31920301
94 -22.07944580 36.36598301
95 -32.03266580 -22.07944580
96 24.50632300 -32.03266580
97 17.55310300 24.50632300
98 -10.40011700 17.55310300
99 -26.84554580 -10.40011700
100 24.69344300 -26.84554580
101 -1.75198581 24.69344300
102 18.29479419 -1.75198581
103 8.34157419 18.29479419
104 9.38835419 8.34157419
> 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/7xsw11290508399.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/8xsw11290508399.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/9xsw11290508399.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/10qjv41290508399.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/11tkcs1290508399.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/12xkay1290508399.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/13bcq71290508399.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/14edpd1290508399.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/15hvni1290508399.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/163w361290508399.tab")
+ }
>
> try(system("convert tmp/11igb1290508399.ps tmp/11igb1290508399.png",intern=TRUE))
character(0)
> try(system("convert tmp/21igb1290508399.ps tmp/21igb1290508399.png",intern=TRUE))
character(0)
> try(system("convert tmp/3uafe1290508399.ps tmp/3uafe1290508399.png",intern=TRUE))
character(0)
> try(system("convert tmp/4uafe1290508399.ps tmp/4uafe1290508399.png",intern=TRUE))
character(0)
> try(system("convert tmp/5uafe1290508399.ps tmp/5uafe1290508399.png",intern=TRUE))
character(0)
> try(system("convert tmp/6mjfy1290508399.ps tmp/6mjfy1290508399.png",intern=TRUE))
character(0)
> try(system("convert tmp/7xsw11290508399.ps tmp/7xsw11290508399.png",intern=TRUE))
character(0)
> try(system("convert tmp/8xsw11290508399.ps tmp/8xsw11290508399.png",intern=TRUE))
character(0)
> try(system("convert tmp/9xsw11290508399.ps tmp/9xsw11290508399.png",intern=TRUE))
character(0)
> try(system("convert tmp/10qjv41290508399.ps tmp/10qjv41290508399.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.480 2.603 6.516