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.
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(147768,0,137507,0,136919,0,136151,0,133001,0,125554,0,119647,0,114158,0,116193,0,152803,0,161761,0,160942,0,149470,0,139208,0,134588,0,130322,0,126611,0,122401,0,117352,0,112135,0,112879,0,148729,0,157230,0,157221,0,146681,0,136524,0,132111,0,125326,1,122716,1,116615,1,113719,1,110737,1,112093,1,143565,1,149946,1,149147,1,134339,1,122683,1,115614,1,116566,1,111272,1,104609,1,101802,1,94542,1,93051,1,124129,1,130374,1,123946,1,114971,1,105531,1,104919,1,104782,0,101281,0,94545,0,93248,0,84031,0,87486,0,115867,0,120327,0,117008,0,108811,0),dim=c(2,61),dimnames=list(c('jonger_dan_25','plan'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('jonger_dan_25','plan'),1:61))
> 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
jonger_dan_25 plan M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 147768 0 1 0 0 0 0 0 0 0 0 0 0 1
2 137507 0 0 1 0 0 0 0 0 0 0 0 0 2
3 136919 0 0 0 1 0 0 0 0 0 0 0 0 3
4 136151 0 0 0 0 1 0 0 0 0 0 0 0 4
5 133001 0 0 0 0 0 1 0 0 0 0 0 0 5
6 125554 0 0 0 0 0 0 1 0 0 0 0 0 6
7 119647 0 0 0 0 0 0 0 1 0 0 0 0 7
8 114158 0 0 0 0 0 0 0 0 1 0 0 0 8
9 116193 0 0 0 0 0 0 0 0 0 1 0 0 9
10 152803 0 0 0 0 0 0 0 0 0 0 1 0 10
11 161761 0 0 0 0 0 0 0 0 0 0 0 1 11
12 160942 0 0 0 0 0 0 0 0 0 0 0 0 12
13 149470 0 1 0 0 0 0 0 0 0 0 0 0 13
14 139208 0 0 1 0 0 0 0 0 0 0 0 0 14
15 134588 0 0 0 1 0 0 0 0 0 0 0 0 15
16 130322 0 0 0 0 1 0 0 0 0 0 0 0 16
17 126611 0 0 0 0 0 1 0 0 0 0 0 0 17
18 122401 0 0 0 0 0 0 1 0 0 0 0 0 18
19 117352 0 0 0 0 0 0 0 1 0 0 0 0 19
20 112135 0 0 0 0 0 0 0 0 1 0 0 0 20
21 112879 0 0 0 0 0 0 0 0 0 1 0 0 21
22 148729 0 0 0 0 0 0 0 0 0 0 1 0 22
23 157230 0 0 0 0 0 0 0 0 0 0 0 1 23
24 157221 0 0 0 0 0 0 0 0 0 0 0 0 24
25 146681 0 1 0 0 0 0 0 0 0 0 0 0 25
26 136524 0 0 1 0 0 0 0 0 0 0 0 0 26
27 132111 0 0 0 1 0 0 0 0 0 0 0 0 27
28 125326 1 0 0 0 1 0 0 0 0 0 0 0 28
29 122716 1 0 0 0 0 1 0 0 0 0 0 0 29
30 116615 1 0 0 0 0 0 1 0 0 0 0 0 30
31 113719 1 0 0 0 0 0 0 1 0 0 0 0 31
32 110737 1 0 0 0 0 0 0 0 1 0 0 0 32
33 112093 1 0 0 0 0 0 0 0 0 1 0 0 33
34 143565 1 0 0 0 0 0 0 0 0 0 1 0 34
35 149946 1 0 0 0 0 0 0 0 0 0 0 1 35
36 149147 1 0 0 0 0 0 0 0 0 0 0 0 36
37 134339 1 1 0 0 0 0 0 0 0 0 0 0 37
38 122683 1 0 1 0 0 0 0 0 0 0 0 0 38
39 115614 1 0 0 1 0 0 0 0 0 0 0 0 39
40 116566 1 0 0 0 1 0 0 0 0 0 0 0 40
41 111272 1 0 0 0 0 1 0 0 0 0 0 0 41
42 104609 1 0 0 0 0 0 1 0 0 0 0 0 42
43 101802 1 0 0 0 0 0 0 1 0 0 0 0 43
44 94542 1 0 0 0 0 0 0 0 1 0 0 0 44
45 93051 1 0 0 0 0 0 0 0 0 1 0 0 45
46 124129 1 0 0 0 0 0 0 0 0 0 1 0 46
47 130374 1 0 0 0 0 0 0 0 0 0 0 1 47
48 123946 1 0 0 0 0 0 0 0 0 0 0 0 48
49 114971 1 1 0 0 0 0 0 0 0 0 0 0 49
50 105531 1 0 1 0 0 0 0 0 0 0 0 0 50
51 104919 1 0 0 1 0 0 0 0 0 0 0 0 51
52 104782 0 0 0 0 1 0 0 0 0 0 0 0 52
53 101281 0 0 0 0 0 1 0 0 0 0 0 0 53
54 94545 0 0 0 0 0 0 1 0 0 0 0 0 54
55 93248 0 0 0 0 0 0 0 1 0 0 0 0 55
56 84031 0 0 0 0 0 0 0 0 1 0 0 0 56
57 87486 0 0 0 0 0 0 0 0 0 1 0 0 57
58 115867 0 0 0 0 0 0 0 0 0 0 1 0 58
59 120327 0 0 0 0 0 0 0 0 0 0 0 1 59
60 117008 0 0 0 0 0 0 0 0 0 0 0 0 60
61 108811 0 1 0 0 0 0 0 0 0 0 0 0 61
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) plan M1 M2 M3 M4
167691.0 2367.8 -11569.6 -20858.1 -23568.9 -25020.1
M5 M6 M7 M8 M9 M10
-27923.7 -33405.5 -36247.2 -41530.6 -39561.2 -6133.4
M11 t
1525.2 -749.6
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-10132.4 -3521.5 490.6 2869.4 9299.4
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 167690.96 2650.70 63.263 < 2e-16 ***
plan 2367.83 1454.93 1.627 0.110328
M1 -11569.57 3091.38 -3.743 0.000496 ***
M2 -20858.12 3248.71 -6.420 6.20e-08 ***
M3 -23568.92 3243.80 -7.266 3.23e-09 ***
M4 -25020.13 3239.41 -7.724 6.60e-10 ***
M5 -27923.74 3235.53 -8.630 2.97e-11 ***
M6 -33405.55 3232.16 -10.335 1.09e-13 ***
M7 -36247.16 3229.31 -11.224 6.76e-15 ***
M8 -41530.57 3226.97 -12.870 < 2e-16 ***
M9 -39561.17 3225.16 -12.266 2.95e-16 ***
M10 -6133.38 3223.86 -1.902 0.063243 .
M11 1525.21 3223.08 0.473 0.638250
t -749.59 40.93 -18.312 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5096 on 47 degrees of freedom
Multiple R-squared: 0.9443, Adjusted R-squared: 0.9289
F-statistic: 61.3 on 13 and 47 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.38908807 0.77817613 0.61091193
[2,] 0.24717386 0.49434772 0.75282614
[3,] 0.17573018 0.35146036 0.82426982
[4,] 0.13526298 0.27052596 0.86473702
[5,] 0.15283610 0.30567221 0.84716390
[6,] 0.13908855 0.27817711 0.86091145
[7,] 0.11904806 0.23809611 0.88095194
[8,] 0.07637570 0.15275140 0.92362430
[9,] 0.06219076 0.12438152 0.93780924
[10,] 0.04740330 0.09480661 0.95259670
[11,] 0.04621682 0.09243364 0.95378318
[12,] 0.04975615 0.09951230 0.95024385
[13,] 0.04695125 0.09390250 0.95304875
[14,] 0.04609113 0.09218227 0.95390887
[15,] 0.11067866 0.22135731 0.88932134
[16,] 0.15206042 0.30412084 0.84793958
[17,] 0.15319715 0.30639429 0.84680285
[18,] 0.11355149 0.22710297 0.88644851
[19,] 0.12841200 0.25682399 0.87158800
[20,] 0.56016071 0.87967858 0.43983929
[21,] 0.75742993 0.48514013 0.24257007
[22,] 0.93448147 0.13103706 0.06551853
[23,] 0.95968634 0.08062731 0.04031366
[24,] 0.97184268 0.05631463 0.02815732
[25,] 0.96498945 0.07002111 0.03501055
[26,] 0.95580726 0.08838547 0.04419274
[27,] 0.90598157 0.18803687 0.09401843
[28,] 0.90186254 0.19627493 0.09813746
> postscript(file="/var/www/html/freestat/rcomp/tmp/1lz1l1229962017.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/2a19g1229962017.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/374cx1229962017.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/4ib8p1229962017.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/57aex1229962017.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 = 61
Frequency = 1
1 2 3 4 5 6
-7603.8039 -7826.6663 -4954.2663 -3521.4663 -3018.2663 -4233.8663
7 8 9 10 11 12
-6549.6663 -6005.6663 -5190.4663 -1258.6663 790.3337 2246.1337
13 14 15 16 17 18
3093.2945 2869.4321 1709.8321 -355.3679 -413.1679 1608.2321
19 20 21 22 23 24
150.4321 966.4321 490.6321 3662.4321 5254.4321 7520.2321
25 26 27 28 29 30
9299.3928 9180.5304 8227.9304 1275.9043 2319.1043 2449.5043
31 32 33 34 35 36
3144.7043 6195.7043 6331.9043 5125.7043 4597.7043 6073.5043
37 38 39 40 41 42
3584.6651 1966.8027 -1641.7973 1511.0027 -129.7973 -561.3973
43 44 45 46 47 48
222.8027 -1004.1973 -3714.9973 -5315.1973 -5979.1973 -10132.3973
49 50 51 52 53 54
-6788.2365 -6190.0989 -3341.6989 1089.9272 1242.1272 737.5272
55 56 57 58 59 60
3031.7272 -152.2728 2082.9272 -2214.2728 -4663.2728 -5707.4728
61
-1585.3120
> postscript(file="/var/www/html/freestat/rcomp/tmp/65i1z1229962017.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 -7603.8039 NA
1 -7826.6663 -7603.8039
2 -4954.2663 -7826.6663
3 -3521.4663 -4954.2663
4 -3018.2663 -3521.4663
5 -4233.8663 -3018.2663
6 -6549.6663 -4233.8663
7 -6005.6663 -6549.6663
8 -5190.4663 -6005.6663
9 -1258.6663 -5190.4663
10 790.3337 -1258.6663
11 2246.1337 790.3337
12 3093.2945 2246.1337
13 2869.4321 3093.2945
14 1709.8321 2869.4321
15 -355.3679 1709.8321
16 -413.1679 -355.3679
17 1608.2321 -413.1679
18 150.4321 1608.2321
19 966.4321 150.4321
20 490.6321 966.4321
21 3662.4321 490.6321
22 5254.4321 3662.4321
23 7520.2321 5254.4321
24 9299.3928 7520.2321
25 9180.5304 9299.3928
26 8227.9304 9180.5304
27 1275.9043 8227.9304
28 2319.1043 1275.9043
29 2449.5043 2319.1043
30 3144.7043 2449.5043
31 6195.7043 3144.7043
32 6331.9043 6195.7043
33 5125.7043 6331.9043
34 4597.7043 5125.7043
35 6073.5043 4597.7043
36 3584.6651 6073.5043
37 1966.8027 3584.6651
38 -1641.7973 1966.8027
39 1511.0027 -1641.7973
40 -129.7973 1511.0027
41 -561.3973 -129.7973
42 222.8027 -561.3973
43 -1004.1973 222.8027
44 -3714.9973 -1004.1973
45 -5315.1973 -3714.9973
46 -5979.1973 -5315.1973
47 -10132.3973 -5979.1973
48 -6788.2365 -10132.3973
49 -6190.0989 -6788.2365
50 -3341.6989 -6190.0989
51 1089.9272 -3341.6989
52 1242.1272 1089.9272
53 737.5272 1242.1272
54 3031.7272 737.5272
55 -152.2728 3031.7272
56 2082.9272 -152.2728
57 -2214.2728 2082.9272
58 -4663.2728 -2214.2728
59 -5707.4728 -4663.2728
60 -1585.3120 -5707.4728
61 NA -1585.3120
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -7826.6663 -7603.8039
[2,] -4954.2663 -7826.6663
[3,] -3521.4663 -4954.2663
[4,] -3018.2663 -3521.4663
[5,] -4233.8663 -3018.2663
[6,] -6549.6663 -4233.8663
[7,] -6005.6663 -6549.6663
[8,] -5190.4663 -6005.6663
[9,] -1258.6663 -5190.4663
[10,] 790.3337 -1258.6663
[11,] 2246.1337 790.3337
[12,] 3093.2945 2246.1337
[13,] 2869.4321 3093.2945
[14,] 1709.8321 2869.4321
[15,] -355.3679 1709.8321
[16,] -413.1679 -355.3679
[17,] 1608.2321 -413.1679
[18,] 150.4321 1608.2321
[19,] 966.4321 150.4321
[20,] 490.6321 966.4321
[21,] 3662.4321 490.6321
[22,] 5254.4321 3662.4321
[23,] 7520.2321 5254.4321
[24,] 9299.3928 7520.2321
[25,] 9180.5304 9299.3928
[26,] 8227.9304 9180.5304
[27,] 1275.9043 8227.9304
[28,] 2319.1043 1275.9043
[29,] 2449.5043 2319.1043
[30,] 3144.7043 2449.5043
[31,] 6195.7043 3144.7043
[32,] 6331.9043 6195.7043
[33,] 5125.7043 6331.9043
[34,] 4597.7043 5125.7043
[35,] 6073.5043 4597.7043
[36,] 3584.6651 6073.5043
[37,] 1966.8027 3584.6651
[38,] -1641.7973 1966.8027
[39,] 1511.0027 -1641.7973
[40,] -129.7973 1511.0027
[41,] -561.3973 -129.7973
[42,] 222.8027 -561.3973
[43,] -1004.1973 222.8027
[44,] -3714.9973 -1004.1973
[45,] -5315.1973 -3714.9973
[46,] -5979.1973 -5315.1973
[47,] -10132.3973 -5979.1973
[48,] -6788.2365 -10132.3973
[49,] -6190.0989 -6788.2365
[50,] -3341.6989 -6190.0989
[51,] 1089.9272 -3341.6989
[52,] 1242.1272 1089.9272
[53,] 737.5272 1242.1272
[54,] 3031.7272 737.5272
[55,] -152.2728 3031.7272
[56,] 2082.9272 -152.2728
[57,] -2214.2728 2082.9272
[58,] -4663.2728 -2214.2728
[59,] -5707.4728 -4663.2728
[60,] -1585.3120 -5707.4728
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -7826.6663 -7603.8039
2 -4954.2663 -7826.6663
3 -3521.4663 -4954.2663
4 -3018.2663 -3521.4663
5 -4233.8663 -3018.2663
6 -6549.6663 -4233.8663
7 -6005.6663 -6549.6663
8 -5190.4663 -6005.6663
9 -1258.6663 -5190.4663
10 790.3337 -1258.6663
11 2246.1337 790.3337
12 3093.2945 2246.1337
13 2869.4321 3093.2945
14 1709.8321 2869.4321
15 -355.3679 1709.8321
16 -413.1679 -355.3679
17 1608.2321 -413.1679
18 150.4321 1608.2321
19 966.4321 150.4321
20 490.6321 966.4321
21 3662.4321 490.6321
22 5254.4321 3662.4321
23 7520.2321 5254.4321
24 9299.3928 7520.2321
25 9180.5304 9299.3928
26 8227.9304 9180.5304
27 1275.9043 8227.9304
28 2319.1043 1275.9043
29 2449.5043 2319.1043
30 3144.7043 2449.5043
31 6195.7043 3144.7043
32 6331.9043 6195.7043
33 5125.7043 6331.9043
34 4597.7043 5125.7043
35 6073.5043 4597.7043
36 3584.6651 6073.5043
37 1966.8027 3584.6651
38 -1641.7973 1966.8027
39 1511.0027 -1641.7973
40 -129.7973 1511.0027
41 -561.3973 -129.7973
42 222.8027 -561.3973
43 -1004.1973 222.8027
44 -3714.9973 -1004.1973
45 -5315.1973 -3714.9973
46 -5979.1973 -5315.1973
47 -10132.3973 -5979.1973
48 -6788.2365 -10132.3973
49 -6190.0989 -6788.2365
50 -3341.6989 -6190.0989
51 1089.9272 -3341.6989
52 1242.1272 1089.9272
53 737.5272 1242.1272
54 3031.7272 737.5272
55 -152.2728 3031.7272
56 2082.9272 -152.2728
57 -2214.2728 2082.9272
58 -4663.2728 -2214.2728
59 -5707.4728 -4663.2728
60 -1585.3120 -5707.4728
> 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/7lc421229962017.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/8323s1229962017.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/96h3n1229962017.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/10s72s1229962017.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/11522p1229962017.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/124se41229962017.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/13e7cw1229962017.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/14dplr1229962017.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/15fsbx1229962017.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/16o1541229962018.tab")
+ }
> system("convert tmp/1lz1l1229962017.ps tmp/1lz1l1229962017.png")
> system("convert tmp/2a19g1229962017.ps tmp/2a19g1229962017.png")
> system("convert tmp/374cx1229962017.ps tmp/374cx1229962017.png")
> system("convert tmp/4ib8p1229962017.ps tmp/4ib8p1229962017.png")
> system("convert tmp/57aex1229962017.ps tmp/57aex1229962017.png")
> system("convert tmp/65i1z1229962017.ps tmp/65i1z1229962017.png")
> system("convert tmp/7lc421229962017.ps tmp/7lc421229962017.png")
> system("convert tmp/8323s1229962017.ps tmp/8323s1229962017.png")
> system("convert tmp/96h3n1229962017.ps tmp/96h3n1229962017.png")
> system("convert tmp/10s72s1229962017.ps tmp/10s72s1229962017.png")
>
>
> proc.time()
user system elapsed
3.582 2.446 4.049