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(9051,0,8823,0,8776,0,8255,0,7969,0,8758,0,8693,0,8271,0,7790,0,7769,0,8170,0,8209,0,9395,0,9260,0,9018,0,8501,0,8500,0,9649,0,9319,0,8830,0,8436,0,8169,0,8269,0,7945,0,9144,0,8770,0,8834,0,7837,0,7792,0,8616,0,8518,0,7940,0,7545,0,7531,0,7665,0,7599,0,8444,0,8549,0,7986,0,7335,0,7287,0,7870,0,7839,0,7327,0,7259,0,6964,0,7271,0,6956,0,7608,0,7692,0,7255,0,6804,0,6655,0,7341,0,7602,0,7086,0,6625,0,6272,0,6576,0,6491,0,7649,0,7400,0,6913,0,6532,0,6486,0,7295,0,7556,0,7088,1,6952,1,6773,1,6917,1,7371,1,8221,1,7953,1,8027,1,7287,1,8076,1,8933,1,9433,1,9479,1,9199,1,9469,1,10015,1,10999,1,13009,1,13699,1,13895,1,13248,1,13973,1,15095,1,15201,1,14823,1,14538,1,14547,1,14407,1),dim=c(2,95),dimnames=list(c('Y','X'),1:95))
> y <- array(NA,dim=c(2,95),dimnames=list(c('Y','X'),1:95))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> #'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
1 9051 0
2 8823 0
3 8776 0
4 8255 0
5 7969 0
6 8758 0
7 8693 0
8 8271 0
9 7790 0
10 7769 0
11 8170 0
12 8209 0
13 9395 0
14 9260 0
15 9018 0
16 8501 0
17 8500 0
18 9649 0
19 9319 0
20 8830 0
21 8436 0
22 8169 0
23 8269 0
24 7945 0
25 9144 0
26 8770 0
27 8834 0
28 7837 0
29 7792 0
30 8616 0
31 8518 0
32 7940 0
33 7545 0
34 7531 0
35 7665 0
36 7599 0
37 8444 0
38 8549 0
39 7986 0
40 7335 0
41 7287 0
42 7870 0
43 7839 0
44 7327 0
45 7259 0
46 6964 0
47 7271 0
48 6956 0
49 7608 0
50 7692 0
51 7255 0
52 6804 0
53 6655 0
54 7341 0
55 7602 0
56 7086 0
57 6625 0
58 6272 0
59 6576 0
60 6491 0
61 7649 0
62 7400 0
63 6913 0
64 6532 0
65 6486 0
66 7295 0
67 7556 0
68 7088 1
69 6952 1
70 6773 1
71 6917 1
72 7371 1
73 8221 1
74 7953 1
75 8027 1
76 7287 1
77 8076 1
78 8933 1
79 9433 1
80 9479 1
81 9199 1
82 9469 1
83 10015 1
84 10999 1
85 13009 1
86 13699 1
87 13895 1
88 13248 1
89 13973 1
90 15095 1
91 15201 1
92 14823 1
93 14538 1
94 14547 1
95 14407 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X
7889 2777
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3892.25 -954.17 -98.67 836.83 4535.75
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7888.7 220.7 35.75 < 2e-16 ***
X 2776.6 406.5 6.83 8.62e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1806 on 93 degrees of freedom
Multiple R-squared: 0.3341, Adjusted R-squared: 0.3269
F-statistic: 46.65 on 1 and 93 DF, p-value: 8.616e-10
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 2.825350e-02 5.650700e-02 0.97174650
[2,] 6.786310e-03 1.357262e-02 0.99321369
[3,] 1.430409e-03 2.860817e-03 0.99856959
[4,] 3.832228e-04 7.664455e-04 0.99961678
[5,] 3.340972e-04 6.681944e-04 0.99966590
[6,] 2.050121e-04 4.100242e-04 0.99979499
[7,] 5.439627e-05 1.087925e-04 0.99994560
[8,] 1.325314e-05 2.650628e-05 0.99998675
[9,] 2.205574e-05 4.411147e-05 0.99997794
[10,] 1.601515e-05 3.203029e-05 0.99998398
[11,] 6.428267e-06 1.285653e-05 0.99999357
[12,] 1.714060e-06 3.428120e-06 0.99999829
[13,] 4.408070e-07 8.816140e-07 0.99999956
[14,] 8.382345e-07 1.676469e-06 0.99999916
[15,] 5.031678e-07 1.006336e-06 0.99999950
[16,] 1.518350e-07 3.036701e-07 0.99999985
[17,] 4.443599e-08 8.887198e-08 0.99999996
[18,] 1.649396e-08 3.298791e-08 0.99999998
[19,] 5.242747e-09 1.048549e-08 0.99999999
[20,] 2.654368e-09 5.308736e-09 1.00000000
[21,] 1.239530e-09 2.479060e-09 1.00000000
[22,] 3.532488e-10 7.064976e-10 1.00000000
[23,] 1.040201e-10 2.080403e-10 1.00000000
[24,] 6.923116e-11 1.384623e-10 1.00000000
[25,] 4.661772e-11 9.323544e-11 1.00000000
[26,] 1.266247e-11 2.532495e-11 1.00000000
[27,] 3.341931e-12 6.683861e-12 1.00000000
[28,] 1.544234e-12 3.088469e-12 1.00000000
[29,] 1.766149e-12 3.532298e-12 1.00000000
[30,] 1.771329e-12 3.542657e-12 1.00000000
[31,] 1.137861e-12 2.275723e-12 1.00000000
[32,] 7.827382e-13 1.565476e-12 1.00000000
[33,] 2.199161e-13 4.398322e-13 1.00000000
[34,] 6.329643e-14 1.265929e-13 1.00000000
[35,] 2.180308e-14 4.360617e-14 1.00000000
[36,] 2.685392e-14 5.370784e-14 1.00000000
[37,] 3.204910e-14 6.409821e-14 1.00000000
[38,] 1.162703e-14 2.325405e-14 1.00000000
[39,] 4.282303e-15 8.564605e-15 1.00000000
[40,] 3.790493e-15 7.580986e-15 1.00000000
[41,] 3.550984e-15 7.101967e-15 1.00000000
[42,] 6.201324e-15 1.240265e-14 1.00000000
[43,] 4.519106e-15 9.038211e-15 1.00000000
[44,] 5.996592e-15 1.199318e-14 1.00000000
[45,] 2.362458e-15 4.724916e-15 1.00000000
[46,] 8.420222e-16 1.684044e-15 1.00000000
[47,] 5.125454e-16 1.025091e-15 1.00000000
[48,] 7.326468e-16 1.465294e-15 1.00000000
[49,] 1.276576e-15 2.553151e-15 1.00000000
[50,] 5.748380e-16 1.149676e-15 1.00000000
[51,] 1.979103e-16 3.958206e-16 1.00000000
[52,] 1.198008e-16 2.396015e-16 1.00000000
[53,] 1.623698e-16 3.247396e-16 1.00000000
[54,] 4.332841e-16 8.665683e-16 1.00000000
[55,] 4.844027e-16 9.688053e-16 1.00000000
[56,] 5.770142e-16 1.154028e-15 1.00000000
[57,] 1.761385e-16 3.522770e-16 1.00000000
[58,] 5.949652e-17 1.189930e-16 1.00000000
[59,] 3.146661e-17 6.293322e-17 1.00000000
[60,] 2.859274e-17 5.718548e-17 1.00000000
[61,] 2.665719e-17 5.331438e-17 1.00000000
[62,] 8.667411e-18 1.733482e-17 1.00000000
[63,] 2.371828e-18 4.743656e-18 1.00000000
[64,] 2.372736e-18 4.745472e-18 1.00000000
[65,] 3.299137e-18 6.598273e-18 1.00000000
[66,] 7.463804e-18 1.492761e-17 1.00000000
[67,] 2.078987e-17 4.157974e-17 1.00000000
[68,] 5.704622e-17 1.140924e-16 1.00000000
[69,] 1.595596e-16 3.191192e-16 1.00000000
[70,] 5.337068e-16 1.067414e-15 1.00000000
[71,] 2.482344e-15 4.964688e-15 1.00000000
[72,] 6.545085e-14 1.309017e-13 1.00000000
[73,] 1.364996e-12 2.729991e-12 1.00000000
[74,] 2.712299e-11 5.424598e-11 1.00000000
[75,] 6.068511e-10 1.213702e-09 1.00000000
[76,] 1.828731e-08 3.657461e-08 0.99999998
[77,] 1.765838e-06 3.531675e-06 0.99999823
[78,] 3.540319e-04 7.080638e-04 0.99964597
[79,] 5.303649e-02 1.060730e-01 0.94696351
[80,] 7.223233e-01 5.553533e-01 0.27767667
[81,] 9.099382e-01 1.801236e-01 0.09006178
[82,] 9.396288e-01 1.207424e-01 0.06037120
[83,] 9.410027e-01 1.179946e-01 0.05899731
[84,] 9.876580e-01 2.468407e-02 0.01234203
[85,] 9.918552e-01 1.628962e-02 0.00814481
[86,] 9.816441e-01 3.671185e-02 0.01835593
> postscript(file="/var/www/html/rcomp/tmp/1twjj1260026500.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/2nem11260026500.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/3vh161260026500.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/49cbw1260026500.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/5rmox1260026500.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 = 95
Frequency = 1
1 2 3 4 5 6
1162.32836 934.32836 887.32836 366.32836 80.32836 869.32836
7 8 9 10 11 12
804.32836 382.32836 -98.67164 -119.67164 281.32836 320.32836
13 14 15 16 17 18
1506.32836 1371.32836 1129.32836 612.32836 611.32836 1760.32836
19 20 21 22 23 24
1430.32836 941.32836 547.32836 280.32836 380.32836 56.32836
25 26 27 28 29 30
1255.32836 881.32836 945.32836 -51.67164 -96.67164 727.32836
31 32 33 34 35 36
629.32836 51.32836 -343.67164 -357.67164 -223.67164 -289.67164
37 38 39 40 41 42
555.32836 660.32836 97.32836 -553.67164 -601.67164 -18.67164
43 44 45 46 47 48
-49.67164 -561.67164 -629.67164 -924.67164 -617.67164 -932.67164
49 50 51 52 53 54
-280.67164 -196.67164 -633.67164 -1084.67164 -1233.67164 -547.67164
55 56 57 58 59 60
-286.67164 -802.67164 -1263.67164 -1616.67164 -1312.67164 -1397.67164
61 62 63 64 65 66
-239.67164 -488.67164 -975.67164 -1356.67164 -1402.67164 -593.67164
67 68 69 70 71 72
-332.67164 -3577.25000 -3713.25000 -3892.25000 -3748.25000 -3294.25000
73 74 75 76 77 78
-2444.25000 -2712.25000 -2638.25000 -3378.25000 -2589.25000 -1732.25000
79 80 81 82 83 84
-1232.25000 -1186.25000 -1466.25000 -1196.25000 -650.25000 333.75000
85 86 87 88 89 90
2343.75000 3033.75000 3229.75000 2582.75000 3307.75000 4429.75000
91 92 93 94 95
4535.75000 4157.75000 3872.75000 3881.75000 3741.75000
> postscript(file="/var/www/html/rcomp/tmp/6pt3l1260026500.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 = 95
Frequency = 1
lag(myerror, k = 1) myerror
0 1162.32836 NA
1 934.32836 1162.32836
2 887.32836 934.32836
3 366.32836 887.32836
4 80.32836 366.32836
5 869.32836 80.32836
6 804.32836 869.32836
7 382.32836 804.32836
8 -98.67164 382.32836
9 -119.67164 -98.67164
10 281.32836 -119.67164
11 320.32836 281.32836
12 1506.32836 320.32836
13 1371.32836 1506.32836
14 1129.32836 1371.32836
15 612.32836 1129.32836
16 611.32836 612.32836
17 1760.32836 611.32836
18 1430.32836 1760.32836
19 941.32836 1430.32836
20 547.32836 941.32836
21 280.32836 547.32836
22 380.32836 280.32836
23 56.32836 380.32836
24 1255.32836 56.32836
25 881.32836 1255.32836
26 945.32836 881.32836
27 -51.67164 945.32836
28 -96.67164 -51.67164
29 727.32836 -96.67164
30 629.32836 727.32836
31 51.32836 629.32836
32 -343.67164 51.32836
33 -357.67164 -343.67164
34 -223.67164 -357.67164
35 -289.67164 -223.67164
36 555.32836 -289.67164
37 660.32836 555.32836
38 97.32836 660.32836
39 -553.67164 97.32836
40 -601.67164 -553.67164
41 -18.67164 -601.67164
42 -49.67164 -18.67164
43 -561.67164 -49.67164
44 -629.67164 -561.67164
45 -924.67164 -629.67164
46 -617.67164 -924.67164
47 -932.67164 -617.67164
48 -280.67164 -932.67164
49 -196.67164 -280.67164
50 -633.67164 -196.67164
51 -1084.67164 -633.67164
52 -1233.67164 -1084.67164
53 -547.67164 -1233.67164
54 -286.67164 -547.67164
55 -802.67164 -286.67164
56 -1263.67164 -802.67164
57 -1616.67164 -1263.67164
58 -1312.67164 -1616.67164
59 -1397.67164 -1312.67164
60 -239.67164 -1397.67164
61 -488.67164 -239.67164
62 -975.67164 -488.67164
63 -1356.67164 -975.67164
64 -1402.67164 -1356.67164
65 -593.67164 -1402.67164
66 -332.67164 -593.67164
67 -3577.25000 -332.67164
68 -3713.25000 -3577.25000
69 -3892.25000 -3713.25000
70 -3748.25000 -3892.25000
71 -3294.25000 -3748.25000
72 -2444.25000 -3294.25000
73 -2712.25000 -2444.25000
74 -2638.25000 -2712.25000
75 -3378.25000 -2638.25000
76 -2589.25000 -3378.25000
77 -1732.25000 -2589.25000
78 -1232.25000 -1732.25000
79 -1186.25000 -1232.25000
80 -1466.25000 -1186.25000
81 -1196.25000 -1466.25000
82 -650.25000 -1196.25000
83 333.75000 -650.25000
84 2343.75000 333.75000
85 3033.75000 2343.75000
86 3229.75000 3033.75000
87 2582.75000 3229.75000
88 3307.75000 2582.75000
89 4429.75000 3307.75000
90 4535.75000 4429.75000
91 4157.75000 4535.75000
92 3872.75000 4157.75000
93 3881.75000 3872.75000
94 3741.75000 3881.75000
95 NA 3741.75000
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 934.32836 1162.32836
[2,] 887.32836 934.32836
[3,] 366.32836 887.32836
[4,] 80.32836 366.32836
[5,] 869.32836 80.32836
[6,] 804.32836 869.32836
[7,] 382.32836 804.32836
[8,] -98.67164 382.32836
[9,] -119.67164 -98.67164
[10,] 281.32836 -119.67164
[11,] 320.32836 281.32836
[12,] 1506.32836 320.32836
[13,] 1371.32836 1506.32836
[14,] 1129.32836 1371.32836
[15,] 612.32836 1129.32836
[16,] 611.32836 612.32836
[17,] 1760.32836 611.32836
[18,] 1430.32836 1760.32836
[19,] 941.32836 1430.32836
[20,] 547.32836 941.32836
[21,] 280.32836 547.32836
[22,] 380.32836 280.32836
[23,] 56.32836 380.32836
[24,] 1255.32836 56.32836
[25,] 881.32836 1255.32836
[26,] 945.32836 881.32836
[27,] -51.67164 945.32836
[28,] -96.67164 -51.67164
[29,] 727.32836 -96.67164
[30,] 629.32836 727.32836
[31,] 51.32836 629.32836
[32,] -343.67164 51.32836
[33,] -357.67164 -343.67164
[34,] -223.67164 -357.67164
[35,] -289.67164 -223.67164
[36,] 555.32836 -289.67164
[37,] 660.32836 555.32836
[38,] 97.32836 660.32836
[39,] -553.67164 97.32836
[40,] -601.67164 -553.67164
[41,] -18.67164 -601.67164
[42,] -49.67164 -18.67164
[43,] -561.67164 -49.67164
[44,] -629.67164 -561.67164
[45,] -924.67164 -629.67164
[46,] -617.67164 -924.67164
[47,] -932.67164 -617.67164
[48,] -280.67164 -932.67164
[49,] -196.67164 -280.67164
[50,] -633.67164 -196.67164
[51,] -1084.67164 -633.67164
[52,] -1233.67164 -1084.67164
[53,] -547.67164 -1233.67164
[54,] -286.67164 -547.67164
[55,] -802.67164 -286.67164
[56,] -1263.67164 -802.67164
[57,] -1616.67164 -1263.67164
[58,] -1312.67164 -1616.67164
[59,] -1397.67164 -1312.67164
[60,] -239.67164 -1397.67164
[61,] -488.67164 -239.67164
[62,] -975.67164 -488.67164
[63,] -1356.67164 -975.67164
[64,] -1402.67164 -1356.67164
[65,] -593.67164 -1402.67164
[66,] -332.67164 -593.67164
[67,] -3577.25000 -332.67164
[68,] -3713.25000 -3577.25000
[69,] -3892.25000 -3713.25000
[70,] -3748.25000 -3892.25000
[71,] -3294.25000 -3748.25000
[72,] -2444.25000 -3294.25000
[73,] -2712.25000 -2444.25000
[74,] -2638.25000 -2712.25000
[75,] -3378.25000 -2638.25000
[76,] -2589.25000 -3378.25000
[77,] -1732.25000 -2589.25000
[78,] -1232.25000 -1732.25000
[79,] -1186.25000 -1232.25000
[80,] -1466.25000 -1186.25000
[81,] -1196.25000 -1466.25000
[82,] -650.25000 -1196.25000
[83,] 333.75000 -650.25000
[84,] 2343.75000 333.75000
[85,] 3033.75000 2343.75000
[86,] 3229.75000 3033.75000
[87,] 2582.75000 3229.75000
[88,] 3307.75000 2582.75000
[89,] 4429.75000 3307.75000
[90,] 4535.75000 4429.75000
[91,] 4157.75000 4535.75000
[92,] 3872.75000 4157.75000
[93,] 3881.75000 3872.75000
[94,] 3741.75000 3881.75000
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 934.32836 1162.32836
2 887.32836 934.32836
3 366.32836 887.32836
4 80.32836 366.32836
5 869.32836 80.32836
6 804.32836 869.32836
7 382.32836 804.32836
8 -98.67164 382.32836
9 -119.67164 -98.67164
10 281.32836 -119.67164
11 320.32836 281.32836
12 1506.32836 320.32836
13 1371.32836 1506.32836
14 1129.32836 1371.32836
15 612.32836 1129.32836
16 611.32836 612.32836
17 1760.32836 611.32836
18 1430.32836 1760.32836
19 941.32836 1430.32836
20 547.32836 941.32836
21 280.32836 547.32836
22 380.32836 280.32836
23 56.32836 380.32836
24 1255.32836 56.32836
25 881.32836 1255.32836
26 945.32836 881.32836
27 -51.67164 945.32836
28 -96.67164 -51.67164
29 727.32836 -96.67164
30 629.32836 727.32836
31 51.32836 629.32836
32 -343.67164 51.32836
33 -357.67164 -343.67164
34 -223.67164 -357.67164
35 -289.67164 -223.67164
36 555.32836 -289.67164
37 660.32836 555.32836
38 97.32836 660.32836
39 -553.67164 97.32836
40 -601.67164 -553.67164
41 -18.67164 -601.67164
42 -49.67164 -18.67164
43 -561.67164 -49.67164
44 -629.67164 -561.67164
45 -924.67164 -629.67164
46 -617.67164 -924.67164
47 -932.67164 -617.67164
48 -280.67164 -932.67164
49 -196.67164 -280.67164
50 -633.67164 -196.67164
51 -1084.67164 -633.67164
52 -1233.67164 -1084.67164
53 -547.67164 -1233.67164
54 -286.67164 -547.67164
55 -802.67164 -286.67164
56 -1263.67164 -802.67164
57 -1616.67164 -1263.67164
58 -1312.67164 -1616.67164
59 -1397.67164 -1312.67164
60 -239.67164 -1397.67164
61 -488.67164 -239.67164
62 -975.67164 -488.67164
63 -1356.67164 -975.67164
64 -1402.67164 -1356.67164
65 -593.67164 -1402.67164
66 -332.67164 -593.67164
67 -3577.25000 -332.67164
68 -3713.25000 -3577.25000
69 -3892.25000 -3713.25000
70 -3748.25000 -3892.25000
71 -3294.25000 -3748.25000
72 -2444.25000 -3294.25000
73 -2712.25000 -2444.25000
74 -2638.25000 -2712.25000
75 -3378.25000 -2638.25000
76 -2589.25000 -3378.25000
77 -1732.25000 -2589.25000
78 -1232.25000 -1732.25000
79 -1186.25000 -1232.25000
80 -1466.25000 -1186.25000
81 -1196.25000 -1466.25000
82 -650.25000 -1196.25000
83 333.75000 -650.25000
84 2343.75000 333.75000
85 3033.75000 2343.75000
86 3229.75000 3033.75000
87 2582.75000 3229.75000
88 3307.75000 2582.75000
89 4429.75000 3307.75000
90 4535.75000 4429.75000
91 4157.75000 4535.75000
92 3872.75000 4157.75000
93 3881.75000 3872.75000
94 3741.75000 3881.75000
> 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/712n31260026500.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/8jasb1260026500.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/9pjt61260026500.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/10tln01260026500.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/11e2na1260026500.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/12e14q1260026500.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/139fy81260026500.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/143thi1260026500.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/15mkrc1260026500.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/16ngwc1260026501.tab")
+ }
>
> system("convert tmp/1twjj1260026500.ps tmp/1twjj1260026500.png")
> system("convert tmp/2nem11260026500.ps tmp/2nem11260026500.png")
> system("convert tmp/3vh161260026500.ps tmp/3vh161260026500.png")
> system("convert tmp/49cbw1260026500.ps tmp/49cbw1260026500.png")
> system("convert tmp/5rmox1260026500.ps tmp/5rmox1260026500.png")
> system("convert tmp/6pt3l1260026500.ps tmp/6pt3l1260026500.png")
> system("convert tmp/712n31260026500.ps tmp/712n31260026500.png")
> system("convert tmp/8jasb1260026500.ps tmp/8jasb1260026500.png")
> system("convert tmp/9pjt61260026500.ps tmp/9pjt61260026500.png")
> system("convert tmp/10tln01260026500.ps tmp/10tln01260026500.png")
>
>
> proc.time()
user system elapsed
2.829 1.629 3.414