graphics.off();           # shuts down all the open graphics plots
objects();                # prints a list of all variables in memory
remove(list=objects());   # sets objects() as a list, and then removes this list
                          # of variables from memory
options(scipen=20);       # sets the decimal places in scientific notation

sink(file = "output.txt", append = FALSE, type = "output", split = TRUE);
# creates a disk file of the terminal/screen output; type can be changed to
# log inputted commands as well.

myData <- read.table("data1.txt", header=TRUE);
# loads a disk file of data series.  the series should be in columns and include
# a header.

y1 <- myData$FOOD/myData$HIBT;
y2 =  myData$CLOTH/myData$HIBT;
# note the assignment operators = or <- can be used interchangably.  However,
# the operator <<- assigns a variable as "global" thus it can be accessed from
# anywhere in the script.
INC   <- log(myData$HIBT);
INCSQ <- INC^2;

X <- cbind(myData$INC, myData$SA, myData$A, myData$AK, myData$SP, myData$INCSQ);
# the "cbind" routine creates a matrix from the columns listed within ().
# note that the SA, A, AK, and SP variables are dummies for different family
# types.
k = nrow(X);
# how many rows are in X?

# SUMMARY STATISTICS ----------------------------------
# Quantiles
print("Quantiles:");
print(summary(X[1:k,1]));
# MAD
print("Mean absolute deviation:");
print(mad(X[1:k,1]));
# MSD
print("Mean squared deviation:");
print( var(X[1:k,1])*((k-1)/k) );
# Variance
print("Variance:");
print(var(X[1:k,1]));
# S.D.
print("Standard deviation:");
print(sd(X[1:k,1]));
# Coefficient of Variation
print("Coefficient of Variation:");
print( sd(X[1:k,1])/mean(X[1:k,1]) );
# Pearson's Coefficient of Variation
print("Pearson's Coefficient of Variation:");
print( 3*(mean(X[1:k,1])-2)/sd(X[1:k,1]) );
# Third moment about the mean
mySkewness <- function(x)
{ 
m3 <- mean((x - mean(x))^3) 
skew <- m3/(sd(x)^3)
skew
}
# this is an example of an R user function.  once defined it can be called at
# any point to retrieve the skewness.
print("Skewness:");
print(mySkewness(X[1:k,1]));

# DO A LINEAR OLS REGRESSION
myLm <- lm(y1 ~ 0 + X);
print(summary(myLm));
myLm <- lm(y2 ~ 0 + X);
print(summary(myLm));
# the "lm" routine regresses y on the columns of X using OLS.  The "0" in front 
# of X instructs the routine to omit the intercept.

res <- residuals(myLm);
# store the residuals from the last regression.
cof <- coefficients(myLm);
# store a vector of estimated coefficients.

fit = myData$INC*cof[1] + myData$SA*cof[2] + myData$A*cof[3] + myData$AK*cof[4] + myData$SP*cof[5] + myData$INCSQ*cof[6];
# create a vector of fitted values.

plot(fit, y2, asp=1);
# plots fit on y2, using a 1:1 aspect ratio

write.table(fit, "fit.txt");
# writes fitted values to a disk file in your working directory.

sink(file = NULL);
# closes the opened disk file opened above with sink