R version 2.9.0 (2009-04-17)
Copyright (C) 2009 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(2360,8.1,2214,7.4,2825,7.3,2355,7.7,2333,8,3016,8,2155,7.7,2172,6.9,2150,6.6,2533,6.9,2058,7.5,2160,7.9,2260,7.7,2498,6.5,2695,6.1,2799,6.4,2947,6.8,2930,7.1,2318,7.3,2540,7.2,2570,7,2669,7,2450,7,2842,7.3,3440,7.5,2678,7.2,2981,7.7,2260,8,2844,7.9,2546,8,2456,8,2295,7.9,2379,7.9,2479,8,2057,8.1,2280,8.1,2351,8.2,2276,8,2548,8.3,2311,8.5,2201,8.6,2725,8.7,2408,8.7,2139,8.5,1898,8.4,2537,8.5,2069,8.7,2063,8.7,2524,8.6,2437,7.9,2189,8.1,2793,8.2,2074,8.5,2622,8.6,2278,8.5,2144,8.3,2427,8.2,2139,8.7,1828,9.3,2072,9.3,1800,8.8,1758,7.4,2246,7.2,1987,7.5,1868,8.3,2514,8.8,2121,8.9),dim=c(2,67),dimnames=list(c('Y','X'),1:67))
> y <- array(NA,dim=c(2,67),dimnames=list(c('Y','X'),1:67))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 2360 8.1 1 0 0 0 0 0 0 0 0 0 0 1
2 2214 7.4 0 1 0 0 0 0 0 0 0 0 0 2
3 2825 7.3 0 0 1 0 0 0 0 0 0 0 0 3
4 2355 7.7 0 0 0 1 0 0 0 0 0 0 0 4
5 2333 8.0 0 0 0 0 1 0 0 0 0 0 0 5
6 3016 8.0 0 0 0 0 0 1 0 0 0 0 0 6
7 2155 7.7 0 0 0 0 0 0 1 0 0 0 0 7
8 2172 6.9 0 0 0 0 0 0 0 1 0 0 0 8
9 2150 6.6 0 0 0 0 0 0 0 0 1 0 0 9
10 2533 6.9 0 0 0 0 0 0 0 0 0 1 0 10
11 2058 7.5 0 0 0 0 0 0 0 0 0 0 1 11
12 2160 7.9 0 0 0 0 0 0 0 0 0 0 0 12
13 2260 7.7 1 0 0 0 0 0 0 0 0 0 0 13
14 2498 6.5 0 1 0 0 0 0 0 0 0 0 0 14
15 2695 6.1 0 0 1 0 0 0 0 0 0 0 0 15
16 2799 6.4 0 0 0 1 0 0 0 0 0 0 0 16
17 2947 6.8 0 0 0 0 1 0 0 0 0 0 0 17
18 2930 7.1 0 0 0 0 0 1 0 0 0 0 0 18
19 2318 7.3 0 0 0 0 0 0 1 0 0 0 0 19
20 2540 7.2 0 0 0 0 0 0 0 1 0 0 0 20
21 2570 7.0 0 0 0 0 0 0 0 0 1 0 0 21
22 2669 7.0 0 0 0 0 0 0 0 0 0 1 0 22
23 2450 7.0 0 0 0 0 0 0 0 0 0 0 1 23
24 2842 7.3 0 0 0 0 0 0 0 0 0 0 0 24
25 3440 7.5 1 0 0 0 0 0 0 0 0 0 0 25
26 2678 7.2 0 1 0 0 0 0 0 0 0 0 0 26
27 2981 7.7 0 0 1 0 0 0 0 0 0 0 0 27
28 2260 8.0 0 0 0 1 0 0 0 0 0 0 0 28
29 2844 7.9 0 0 0 0 1 0 0 0 0 0 0 29
30 2546 8.0 0 0 0 0 0 1 0 0 0 0 0 30
31 2456 8.0 0 0 0 0 0 0 1 0 0 0 0 31
32 2295 7.9 0 0 0 0 0 0 0 1 0 0 0 32
33 2379 7.9 0 0 0 0 0 0 0 0 1 0 0 33
34 2479 8.0 0 0 0 0 0 0 0 0 0 1 0 34
35 2057 8.1 0 0 0 0 0 0 0 0 0 0 1 35
36 2280 8.1 0 0 0 0 0 0 0 0 0 0 0 36
37 2351 8.2 1 0 0 0 0 0 0 0 0 0 0 37
38 2276 8.0 0 1 0 0 0 0 0 0 0 0 0 38
39 2548 8.3 0 0 1 0 0 0 0 0 0 0 0 39
40 2311 8.5 0 0 0 1 0 0 0 0 0 0 0 40
41 2201 8.6 0 0 0 0 1 0 0 0 0 0 0 41
42 2725 8.7 0 0 0 0 0 1 0 0 0 0 0 42
43 2408 8.7 0 0 0 0 0 0 1 0 0 0 0 43
44 2139 8.5 0 0 0 0 0 0 0 1 0 0 0 44
45 1898 8.4 0 0 0 0 0 0 0 0 1 0 0 45
46 2537 8.5 0 0 0 0 0 0 0 0 0 1 0 46
47 2069 8.7 0 0 0 0 0 0 0 0 0 0 1 47
48 2063 8.7 0 0 0 0 0 0 0 0 0 0 0 48
49 2524 8.6 1 0 0 0 0 0 0 0 0 0 0 49
50 2437 7.9 0 1 0 0 0 0 0 0 0 0 0 50
51 2189 8.1 0 0 1 0 0 0 0 0 0 0 0 51
52 2793 8.2 0 0 0 1 0 0 0 0 0 0 0 52
53 2074 8.5 0 0 0 0 1 0 0 0 0 0 0 53
54 2622 8.6 0 0 0 0 0 1 0 0 0 0 0 54
55 2278 8.5 0 0 0 0 0 0 1 0 0 0 0 55
56 2144 8.3 0 0 0 0 0 0 0 1 0 0 0 56
57 2427 8.2 0 0 0 0 0 0 0 0 1 0 0 57
58 2139 8.7 0 0 0 0 0 0 0 0 0 1 0 58
59 1828 9.3 0 0 0 0 0 0 0 0 0 0 1 59
60 2072 9.3 0 0 0 0 0 0 0 0 0 0 0 60
61 1800 8.8 1 0 0 0 0 0 0 0 0 0 0 61
62 1758 7.4 0 1 0 0 0 0 0 0 0 0 0 62
63 2246 7.2 0 0 1 0 0 0 0 0 0 0 0 63
64 1987 7.5 0 0 0 1 0 0 0 0 0 0 0 64
65 1868 8.3 0 0 0 0 1 0 0 0 0 0 0 65
66 2514 8.8 0 0 0 0 0 1 0 0 0 0 0 66
67 2121 8.9 0 0 0 0 0 0 1 0 0 0 0 67
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
3519.2088 -133.2486 139.0016 -102.8466 178.0707 54.1918
M5 M6 M7 M8 M9 M10
58.2546 434.1051 -0.5275 -107.0438 -95.1437 121.8609
M11 t
-213.4097 -3.7549
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-456.58 -174.60 17.30 134.99 875.03
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3519.2088 557.4624 6.313 5.73e-08 ***
X -133.2486 72.7907 -1.831 0.0728 .
M1 139.0016 165.7534 0.839 0.4055
M2 -102.8466 174.8244 -0.588 0.5588
M3 178.0707 174.1522 1.023 0.3112
M4 54.1918 169.4656 0.320 0.7504
M5 58.2546 166.3081 0.350 0.7275
M6 434.1051 165.5833 2.622 0.0114 *
M7 -0.5275 165.6945 -0.003 0.9975
M8 -107.0438 175.5426 -0.610 0.5446
M9 -95.1437 177.8747 -0.535 0.5950
M10 121.8609 175.2767 0.695 0.4899
M11 -213.4097 173.1034 -1.233 0.2231
t -3.7549 2.4239 -1.549 0.1273
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 273.4 on 53 degrees of freedom
Multiple R-squared: 0.4403, Adjusted R-squared: 0.303
F-statistic: 3.207 on 13 and 53 DF, p-value: 0.001322
> 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.4968054 0.9936108 0.50319460
[2,] 0.3969383 0.7938766 0.60306172
[3,] 0.3474070 0.6948139 0.65259305
[4,] 0.3568984 0.7137969 0.64310156
[5,] 0.2793958 0.5587916 0.72060418
[6,] 0.1967644 0.3935288 0.80323560
[7,] 0.1404703 0.2809406 0.85952970
[8,] 0.1858067 0.3716135 0.81419327
[9,] 0.6452322 0.7095355 0.35476775
[10,] 0.6223880 0.7552240 0.37761198
[11,] 0.6335943 0.7328113 0.36640566
[12,] 0.8102119 0.3795762 0.18978812
[13,] 0.8635009 0.2729983 0.13649913
[14,] 0.9376580 0.1246840 0.06234200
[15,] 0.9055786 0.1888427 0.09442137
[16,] 0.8662062 0.2675876 0.13379381
[17,] 0.8098072 0.3803856 0.19019280
[18,] 0.7492300 0.5015401 0.25077004
[19,] 0.6964886 0.6070229 0.30351144
[20,] 0.6511149 0.6977703 0.34888515
[21,] 0.6567967 0.6864066 0.34320328
[22,] 0.5787054 0.8425893 0.42129465
[23,] 0.4979879 0.9959757 0.50201213
[24,] 0.4838155 0.9676310 0.51618449
[25,] 0.4468406 0.8936812 0.55315938
[26,] 0.3612132 0.7224265 0.63878676
[27,] 0.2889544 0.5779088 0.71104562
[28,] 0.2308115 0.4616230 0.76918849
[29,] 0.5726648 0.8546705 0.42733523
[30,] 0.4692542 0.9385083 0.53074583
[31,] 0.3485387 0.6970773 0.65146133
[32,] 0.2954716 0.5909432 0.70452840
[33,] 0.3242969 0.6485937 0.67570314
[34,] 0.3003610 0.6007220 0.69963901
> postscript(file="/var/www/html/rcomp/tmp/18cwh1258742973.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/2j9wm1258742973.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/3jsje1258742973.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/480td1258742973.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/5oujf1258742973.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 = 67
Frequency = 1
1 2 3 4 5 6
-215.1418398 -208.8127424 111.6999655 -177.3668860 -159.7002193 151.2042014
7 8 9 10 11 12
-311.3829059 -290.7106688 -360.8304471 -151.1055857 -207.1310014 -261.4863617
13 14 15 16 17 18
-323.3828552 0.3219351 -133.1399413 138.4683458 339.4598739 -9.6611211
19 20 21 22 23 24
-156.6239213 162.3223457 157.5274289 43.2777060 163.3031218 385.6229001
25 26 27 28 29 30
875.0258523 318.6543953 411.1162717 -142.2754411 428.0917799 -228.6789380
31 32 33 34 35 36
119.7085389 55.6548060 131.5096120 31.5847506 -38.0649723 -24.7197783
37 38 39 40 41 42
-75.6416875 68.3117170 103.1238706 20.4072963 -76.5757599 88.6535223
43 44 45 46 47 48
210.0409992 24.6624048 -237.8076506 201.2674880 98.9426265 -116.7121795
49 50 51 52 53 54
195.7161885 261.0452859 -237.4674220 507.4911423 -171.8421910 17.3870911
55 56 57 58 59 60
98.4497066 48.0711123 309.6010568 -125.0243589 -17.0497746 17.2954194
61 62 63 64 65 66
-456.5756584 -439.5205909 -255.3327445 -346.7244574 -359.4334836 -18.9047557
67
39.8075826
> postscript(file="/var/www/html/rcomp/tmp/6mqeu1258742973.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 = 67
Frequency = 1
lag(myerror, k = 1) myerror
0 -215.1418398 NA
1 -208.8127424 -215.1418398
2 111.6999655 -208.8127424
3 -177.3668860 111.6999655
4 -159.7002193 -177.3668860
5 151.2042014 -159.7002193
6 -311.3829059 151.2042014
7 -290.7106688 -311.3829059
8 -360.8304471 -290.7106688
9 -151.1055857 -360.8304471
10 -207.1310014 -151.1055857
11 -261.4863617 -207.1310014
12 -323.3828552 -261.4863617
13 0.3219351 -323.3828552
14 -133.1399413 0.3219351
15 138.4683458 -133.1399413
16 339.4598739 138.4683458
17 -9.6611211 339.4598739
18 -156.6239213 -9.6611211
19 162.3223457 -156.6239213
20 157.5274289 162.3223457
21 43.2777060 157.5274289
22 163.3031218 43.2777060
23 385.6229001 163.3031218
24 875.0258523 385.6229001
25 318.6543953 875.0258523
26 411.1162717 318.6543953
27 -142.2754411 411.1162717
28 428.0917799 -142.2754411
29 -228.6789380 428.0917799
30 119.7085389 -228.6789380
31 55.6548060 119.7085389
32 131.5096120 55.6548060
33 31.5847506 131.5096120
34 -38.0649723 31.5847506
35 -24.7197783 -38.0649723
36 -75.6416875 -24.7197783
37 68.3117170 -75.6416875
38 103.1238706 68.3117170
39 20.4072963 103.1238706
40 -76.5757599 20.4072963
41 88.6535223 -76.5757599
42 210.0409992 88.6535223
43 24.6624048 210.0409992
44 -237.8076506 24.6624048
45 201.2674880 -237.8076506
46 98.9426265 201.2674880
47 -116.7121795 98.9426265
48 195.7161885 -116.7121795
49 261.0452859 195.7161885
50 -237.4674220 261.0452859
51 507.4911423 -237.4674220
52 -171.8421910 507.4911423
53 17.3870911 -171.8421910
54 98.4497066 17.3870911
55 48.0711123 98.4497066
56 309.6010568 48.0711123
57 -125.0243589 309.6010568
58 -17.0497746 -125.0243589
59 17.2954194 -17.0497746
60 -456.5756584 17.2954194
61 -439.5205909 -456.5756584
62 -255.3327445 -439.5205909
63 -346.7244574 -255.3327445
64 -359.4334836 -346.7244574
65 -18.9047557 -359.4334836
66 39.8075826 -18.9047557
67 NA 39.8075826
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -208.8127424 -215.1418398
[2,] 111.6999655 -208.8127424
[3,] -177.3668860 111.6999655
[4,] -159.7002193 -177.3668860
[5,] 151.2042014 -159.7002193
[6,] -311.3829059 151.2042014
[7,] -290.7106688 -311.3829059
[8,] -360.8304471 -290.7106688
[9,] -151.1055857 -360.8304471
[10,] -207.1310014 -151.1055857
[11,] -261.4863617 -207.1310014
[12,] -323.3828552 -261.4863617
[13,] 0.3219351 -323.3828552
[14,] -133.1399413 0.3219351
[15,] 138.4683458 -133.1399413
[16,] 339.4598739 138.4683458
[17,] -9.6611211 339.4598739
[18,] -156.6239213 -9.6611211
[19,] 162.3223457 -156.6239213
[20,] 157.5274289 162.3223457
[21,] 43.2777060 157.5274289
[22,] 163.3031218 43.2777060
[23,] 385.6229001 163.3031218
[24,] 875.0258523 385.6229001
[25,] 318.6543953 875.0258523
[26,] 411.1162717 318.6543953
[27,] -142.2754411 411.1162717
[28,] 428.0917799 -142.2754411
[29,] -228.6789380 428.0917799
[30,] 119.7085389 -228.6789380
[31,] 55.6548060 119.7085389
[32,] 131.5096120 55.6548060
[33,] 31.5847506 131.5096120
[34,] -38.0649723 31.5847506
[35,] -24.7197783 -38.0649723
[36,] -75.6416875 -24.7197783
[37,] 68.3117170 -75.6416875
[38,] 103.1238706 68.3117170
[39,] 20.4072963 103.1238706
[40,] -76.5757599 20.4072963
[41,] 88.6535223 -76.5757599
[42,] 210.0409992 88.6535223
[43,] 24.6624048 210.0409992
[44,] -237.8076506 24.6624048
[45,] 201.2674880 -237.8076506
[46,] 98.9426265 201.2674880
[47,] -116.7121795 98.9426265
[48,] 195.7161885 -116.7121795
[49,] 261.0452859 195.7161885
[50,] -237.4674220 261.0452859
[51,] 507.4911423 -237.4674220
[52,] -171.8421910 507.4911423
[53,] 17.3870911 -171.8421910
[54,] 98.4497066 17.3870911
[55,] 48.0711123 98.4497066
[56,] 309.6010568 48.0711123
[57,] -125.0243589 309.6010568
[58,] -17.0497746 -125.0243589
[59,] 17.2954194 -17.0497746
[60,] -456.5756584 17.2954194
[61,] -439.5205909 -456.5756584
[62,] -255.3327445 -439.5205909
[63,] -346.7244574 -255.3327445
[64,] -359.4334836 -346.7244574
[65,] -18.9047557 -359.4334836
[66,] 39.8075826 -18.9047557
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -208.8127424 -215.1418398
2 111.6999655 -208.8127424
3 -177.3668860 111.6999655
4 -159.7002193 -177.3668860
5 151.2042014 -159.7002193
6 -311.3829059 151.2042014
7 -290.7106688 -311.3829059
8 -360.8304471 -290.7106688
9 -151.1055857 -360.8304471
10 -207.1310014 -151.1055857
11 -261.4863617 -207.1310014
12 -323.3828552 -261.4863617
13 0.3219351 -323.3828552
14 -133.1399413 0.3219351
15 138.4683458 -133.1399413
16 339.4598739 138.4683458
17 -9.6611211 339.4598739
18 -156.6239213 -9.6611211
19 162.3223457 -156.6239213
20 157.5274289 162.3223457
21 43.2777060 157.5274289
22 163.3031218 43.2777060
23 385.6229001 163.3031218
24 875.0258523 385.6229001
25 318.6543953 875.0258523
26 411.1162717 318.6543953
27 -142.2754411 411.1162717
28 428.0917799 -142.2754411
29 -228.6789380 428.0917799
30 119.7085389 -228.6789380
31 55.6548060 119.7085389
32 131.5096120 55.6548060
33 31.5847506 131.5096120
34 -38.0649723 31.5847506
35 -24.7197783 -38.0649723
36 -75.6416875 -24.7197783
37 68.3117170 -75.6416875
38 103.1238706 68.3117170
39 20.4072963 103.1238706
40 -76.5757599 20.4072963
41 88.6535223 -76.5757599
42 210.0409992 88.6535223
43 24.6624048 210.0409992
44 -237.8076506 24.6624048
45 201.2674880 -237.8076506
46 98.9426265 201.2674880
47 -116.7121795 98.9426265
48 195.7161885 -116.7121795
49 261.0452859 195.7161885
50 -237.4674220 261.0452859
51 507.4911423 -237.4674220
52 -171.8421910 507.4911423
53 17.3870911 -171.8421910
54 98.4497066 17.3870911
55 48.0711123 98.4497066
56 309.6010568 48.0711123
57 -125.0243589 309.6010568
58 -17.0497746 -125.0243589
59 17.2954194 -17.0497746
60 -456.5756584 17.2954194
61 -439.5205909 -456.5756584
62 -255.3327445 -439.5205909
63 -346.7244574 -255.3327445
64 -359.4334836 -346.7244574
65 -18.9047557 -359.4334836
66 39.8075826 -18.9047557
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/7n0t11258742973.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/839e81258742973.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/92mi71258742973.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/rcomp/tmp/108kt21258742973.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/111pbd1258742973.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/12ea1w1258742973.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/137nju1258742973.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/14gu3b1258742973.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/1540jh1258742973.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/168x3k1258742973.tab")
+ }
>
> system("convert tmp/18cwh1258742973.ps tmp/18cwh1258742973.png")
> system("convert tmp/2j9wm1258742973.ps tmp/2j9wm1258742973.png")
> system("convert tmp/3jsje1258742973.ps tmp/3jsje1258742973.png")
> system("convert tmp/480td1258742973.ps tmp/480td1258742973.png")
> system("convert tmp/5oujf1258742973.ps tmp/5oujf1258742973.png")
> system("convert tmp/6mqeu1258742973.ps tmp/6mqeu1258742973.png")
> system("convert tmp/7n0t11258742973.ps tmp/7n0t11258742973.png")
> system("convert tmp/839e81258742973.ps tmp/839e81258742973.png")
> system("convert tmp/92mi71258742973.ps tmp/92mi71258742973.png")
> system("convert tmp/108kt21258742973.ps tmp/108kt21258742973.png")
>
>
> proc.time()
user system elapsed
2.542 1.605 2.936