ECON457 lab02 R For Economics
Jon
September 23, 2016
A series of video of R tutorial.
R for ECON
R for ECON345
library(foreign)
wage1 <- read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/wage1.dta")
# OLS regression:
ols.model <- lm(log(wage) ~ educ+exper+tenure, data=wage1) # Different functions to work with the model
# Print a summary of the results
summary(ols.model)##
## Call:
## lm(formula = log(wage) ~ educ + exper + tenure, data = wage1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.05802 -0.29645 -0.03265 0.28788 1.42809
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.284360 0.104190 2.729 0.00656 **
## educ 0.092029 0.007330 12.555 < 2e-16 ***
## exper 0.004121 0.001723 2.391 0.01714 *
## tenure 0.022067 0.003094 7.133 3.29e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4409 on 522 degrees of freedom
## Multiple R-squared: 0.316, Adjusted R-squared: 0.3121
## F-statistic: 80.39 on 3 and 522 DF, p-value: < 2.2e-16# Produce confidence intervals
confint(ols.model)## 2.5 % 97.5 %
## (Intercept) 0.0796755842 0.48904353
## educ 0.0776292137 0.10642876
## exper 0.0007356983 0.00750652
## tenure 0.0159896850 0.02814475# Get a vector of the coefficients
coef(ols.model)## (Intercept) educ exper tenure
## 0.284359555 0.092028987 0.004121109 0.022067217R for ECON345/365/366
Matrix transposition
M <- cbind(c(1,0,1),c(0,1,2),c(0,0,1))
t(M)## [,1] [,2] [,3]
## [1,] 1 0 1
## [2,] 0 1 2
## [3,] 0 0 1Matrix inversion
solve(M)## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] -1 -2 1solve(M)%*%M## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1R for ECON208
Solving a linear equation
m=matrix(nrow=2,ncol=2,c(1,-.8,1,.2))
l=matrix(c(1.0+25.0/18,25.0/18.0))
(k=solve(m,l))## [,1]
## [1,] -0.9111111
## [2,] 3.3000000m%*%k #checking the answer## [,1]
## [1,] 2.388889
## [2,] 1.388889Eigenvalues and eigenvectors
eigen(m)## $values
## [1] 0.6+0.8i 0.6-0.8i
##
## $vectors
## [,1] [,2]
## [1,] 0.7453560+0.0000000i 0.7453560+0.0000000i
## [2,] -0.2981424+0.5962848i -0.2981424-0.5962848iR for ECON350/351
Symbolic calculations
D(expression(x^n),"x")## x^(n - 1) * nD(expression(exp(a*x)),"x")## exp(a * x) * aCompute derivatives of simple expressions, symbolically and algorithmically.
trig.exp <- expression(sin(cos(x + y^2)))
( D.sc <- D(trig.exp, "x") )## -(cos(cos(x + y^2)) * sin(x + y^2))all.equal(D(trig.exp[[1]], "x"), D.sc)## [1] TRUE## formula argument :
(MU1 <- deriv(~ x1^0.25*x2^0.25, "x1" )) ; ## expression({
## .expr2 <- x2^0.25
## .value <- x1^0.25 * .expr2
## .grad <- array(0, c(length(.value), 1L), list(NULL, c("x1")))
## .grad[, "x1"] <- 0.25 * x1^-0.75 * .expr2
## attr(.value, "gradient") <- .grad
## .value
## })(MU2<- deriv(~ x1^0.25*x2^0.25, "x1" ))## expression({
## .expr2 <- x2^0.25
## .value <- x1^0.25 * .expr2
## .grad <- array(0, c(length(.value), 1L), list(NULL, c("x1")))
## .grad[, "x1"] <- 0.25 * x1^-0.75 * .expr2
## attr(.value, "gradient") <- .grad
## .value
## })mode(MU1)## [1] "expression"x1<- 2
x2 <- 2
eval(MU1)## [1] 1.414214
## attr(,"gradient")
## x1
## [1,] 0.1767767Unconstrained optimization
f <- function(x){(x[1] - 5)^2 + (x[2] - 6)^2}
initial_x <- c(10, 11)
x_optimal <- optim(initial_x, f, method="CG")
x_min <- x_optimal$par
x_min## [1] 5 6ECON457
Equality constraint optimization
Inequality constraint optimization
Linear Programming
Maximize expected return: f(x1, x2, x3) = x15% + x24% + x3*6%
Subjected to constraints: 10% < x1, x2, x3 < 100% x1 + x2 + x3 = 1 x3 < x1 + x2 x1 < 2 * x2
#install.packages("lpSolve")
#install.packages("lpSolveAPI")
library(lpSolve)
library(lpSolveAPI)
# Set the number of vars
model <- make.lp(0, 3)
# Define the object function: for Minimize, use -ve
set.objfn(model, c(-0.05, -0.04, -0.06))
# Add the constraints
add.constraint(model, c(1, 1, 1), "=", 1)
add.constraint(model, c(1, 1, -1), ">", 0)
add.constraint(model, c(1, -2, 0), "<", 0)
# Set the upper and lower bounds
set.bounds(model, lower=c(0.1, 0.1, 0.1), upper=c(1, 1, 1))
# Compute the optimized model
solve(model)## [1] 0# Get the value of the optimized parameters
get.variables(model)## [1] 0.3333333 0.1666667 0.5000000# Get the value of the objective function
get.objective(model)## [1] -0.05333333# Get the value of the constraint
get.constraints(model)## [1] 1 0 0Quadratic Programming
To formulate the problem:
Variable: x1 = % investment in A, x2 = % in B, x3 = % in C
Minimize variance: x
Subjected to constraint: x1 + x2 + x3 == 1 X15% + x24% + x3*6% >= 5.2% 0% < x1, x2, x3 < 100%
#install.packages("quadprog")
library(quadprog)
mu_return_vector <- c(0.05, 0.04, 0.06)
sigma <- matrix(c(0.01, 0.002, 0.005,
0.002, 0.008, 0.006,
0.005, 0.006, 0.012),
nrow=3, ncol=3)
D.Matrix <- 2*sigma
d.Vector <- rep(0, 3)
A.Equality <- matrix(c(1,1,1), ncol=1)
A.Matrix <- cbind(A.Equality, mu_return_vector,
diag(3))
b.Vector <- c(1, 0.052, rep(0, 3))
out <- solve.QP(Dmat=D.Matrix, dvec=d.Vector,
Amat=A.Matrix, bvec=b.Vector,
meq=1)
out$solution## [1] 0.4 0.2 0.4out$value## [1] 0.00672Ryacas: Symbolic calculation
only works for R 3.3.2
# # options(repos = c(CRAN = "http://cran.rstudio.com"))
# # install.packages('Ryacas')
# library(Ryacas)
#
# #package?Ryacas
# print(yacas(expression(integrate(1/x, x))))
# print(yacas("Integrate(x) 1/x"))
# x <- Sym("x"); Integrate(1/x, x)
# acos(Sym("1/2"))
# my_func <- function(x) {
# return(x/(x^2 + 3))
# }
# my_deriv <- yacas(deriv(my_func(x), x))R for ECON435/468
http://www.quantmod.com/ https://www.youtube.com/watch?v=61_F2fcvrsw
#install.packages("quantmod")
library(quantmod)## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: TTR## Version 0.4-0 included new data defaults. See ?getSymbols.getSymbols(c("AMZN"))## As of 0.4-0, 'getSymbols' uses env=parent.frame() and
## auto.assign=TRUE by default.
##
## This behavior will be phased out in 0.5-0 when the call will
## default to use auto.assign=FALSE. getOption("getSymbols.env") and
## getOptions("getSymbols.auto.assign") are now checked for alternate defaults
##
## This message is shown once per session and may be disabled by setting
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for more details.## [1] "AMZN"barChart(AMZN,theme='white.mono',bar.type='hlc')
getSymbols(c("^GSPC"))## [1] "GSPC"chartSeries(GSPC, subset='last 3 months')
addBBands(n = 20, sd = 2, ma = "SMA", draw = 'bands', on = -1)
See what packages we used.
sessionInfo()## R version 3.3.1 (2016-06-21)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 7 x64 (build 7601) Service Pack 1
##
## locale:
## [1] LC_COLLATE=English_Canada.1252 LC_CTYPE=English_Canada.1252
## [3] LC_MONETARY=English_Canada.1252 LC_NUMERIC=C
## [5] LC_TIME=English_Canada.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] quantmod_0.4-6 TTR_0.23-1 xts_0.9-7
## [4] zoo_1.7-13 quadprog_1.5-5 lpSolveAPI_5.5.2.0-17
## [7] lpSolve_5.6.13 foreign_0.8-67
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.7 lattice_0.20-34 digest_0.6.10 grid_3.3.1
## [5] magrittr_1.5 evaluate_0.9 stringi_1.1.1 rmarkdown_1.0
## [9] tools_3.3.1 stringr_1.1.0 yaml_2.1.13 htmltools_0.3.5
## [13] knitr_1.14