# lecturenotes10.r

rm(list = ls(all = TRUE))

setwd("C:\\Users\\Kevin\\Documents\\CCNY\\data for classes\\CPS2010")

data1 <-read.spss("cps_m2010_work3.sav", use.value.labels=T, to.data.frame=T)
str(data1)

data1$WEXP <- as.numeric(data1$WEXP)

# restrict to prime-age 25-55 and IncomeWage > 0 and FT worker
restrict1 <- ((data1$A_AGE >= 25) & (data1$A_AGE <= 55) & (data1$WSAL_VAL > 0) & ((data1$WEXP == 1) | (data1$WEXP  == 2)  | (data1$WEXP  == 3) ))


data2 <- subset(data1, restrict1)
str(data2)

restrict2 <- as.logical(data2$ed_coll)
data3 <- subset(data2, restrict2)

Age25 <- data3$WSAL_VAL[data3$A_AGE == 25]
Age30 <- data3$WSAL_VAL[data3$A_AGE == 30]
Age35 <- data3$WSAL_VAL[data3$A_AGE == 35]
Age40 <- data3$WSAL_VAL[data3$A_AGE == 40]
Age45 <- data3$WSAL_VAL[data3$A_AGE == 45]
Age50 <- data3$WSAL_VAL[data3$A_AGE == 50]

plot(density(Age25, bw = "sj"))

# can use nonparametric models - this is the "np" package in R
# follow along with that "vignette" except with our dataset

NN <- length(data3$WSAL_VAL)
restrict3 <- as.logical(round(runif(NN,min=0,max=0.75)))
data4 <- subset(data3, restrict3)


library(np)

model_parametric1 <- lm(WSAL_VAL ~ A_AGE + I(A_AGE^2), data = data3)
summary(model_parametric1)

# note that this is rather computationally intensive!
model_nonparametric1 <- npreg(WSAL_VAL ~ A_AGE, regtype = "ll", bwmethod = "cv.aic", gradients = TRUE, data = data3)
summary(model_nonparametric1)

npsigtest(model_nonparametric1)


plot(data3$A_AGE,data3$WSAL_VAL, xlab = "age", ylab = "wage", cex=.1)
lines(data3$A_AGE, fitted(model_parametric1), lty = 2, col = "red")
lines(data3$A_AGE, fitted(model_nonparametric1), lty = 1, col = "blue")

# ------------------
# time series
# ------------------
library(lattice)
library(latticeExtra)
rm(list = ls(all = TRUE))

setwd("C:\\Users\\Kevin\\Documents\\CCNY\\data for classes")
dat_empl_NSA = read.csv("empl_data.csv")
# Year,Month,CivLabForceNSA,EmplNSA,UnemplNSA
# 1948,1,58690,56339,2351
# 1948,2,59247,56440,2807
dat_empl_SA = read.csv("empl_dataSA.csv")
# Year,Month,CivLabForSA,Empl SA,Unempl SA
# 1948,1,60095,58061,2034
# 1948,2,60524,58196,2328

CivLabForceNSA <- ts(dat_empl_NSA$CivLabForceNSA, frequency = 12, start = c(1948,1))
CivLabForce <- ts(dat_empl_SA$CivLabForSA, frequency = 12, start = c(1948,1))

asTheEconomist(xyplot(CivLabForceNSA, aspect=0.25, ylab="Civilian Labor Force NSA"))
asTheEconomist(xyplot(CivLabForce, aspect=0.25, ylab="Civilian Labor Force"))

ar1 <- arima(CivLabForceNSA,order=c(1,0,0))
tsdiag(ar1)
yr_ahead1 <- predict(ar1, n.ahead=12)

ma1 <- arima(CivLabForceNSA,order=c(0,0,1))
tsdiag(ma1)

ar2 <- arima(CivLabForce,order=c(1,0,0))
tsdiag(ar2)
ma2 <- arima(CivLabForce,order=c(0,0,1))
tsdiag(ma2)

CLF_d1 <- diff(CivLabForce)
asTheEconomist(xyplot(CLF_d1, aspect=0.25, ylab="diff Civilian Labor Force"))
ar3 <- arima(CLF_d1, order=c(1,0,0))
tsdiag(ar3)

CLFNSA_d1 <- diff(CivLabForceNSA)
asTheEconomist(xyplot(CLFNSA_d1, aspect=0.25, ylab="diff Civilian Labor Force"))
ar4 <- arima(CLFNSA_d1, order=c(1,0,0))
tsdiag(ar4)