## Data input & Data output 
## Loops
1. for loop
2. while loop
3. Ifelse & break
## Power Example

#==========================================================================================
#=========================================================================================
## Data input & Data output
setwd("C:/Users/user/Documents/My Dropbox/Work-R lecture/2012.09.24") # Change Direction
Data <-read.table("Unemployment_Data.txt", header=T)                  # Input data
Output <-Data[,2:3]                           
write.table(Output, file = "data_out.csv", sep = ",", col.names = NA) # print output to a file


## Loops
# Loops Syntax
for(variable in sequence) { 
	statements 
}

# Loops Example 1
# Create a Fibonacci Sequence (1, 1, 2, 3, 5, 8, ...)
num <-10
Fibonacci <-numeric(num)
Fibonacci[1] <-1
Fibonacci[2] <-1

for (i in 3:num){
	Fibonacci[i] <-Fibonacci[i-1]+Fibonacci[i-2]
}
Fibonacci


# Loops Example 2
for (i in 1:5){
	for (j in 1:3){
		print (c(i,j))
	}
}


# While Loop Syntax
while(condition) {
	statements
}

# While Loop Example
z <- 0       
while(z<5) { 
	z <- z+2
	print(z)  
}



# Ifelse & break
# Convergence
n <-1
x <-1/n
repeat{
if(x<0.001){
	break
}else{
	n <-n+1    
	x <-1/n
}
print(n)
print(x)
}



# function Syntax
function(arg1, arg2, ... ){
	statements
	return(object)
}

## Power Example
# Refer to page 15 in the lecture note (Getting Started With R)
# DGP: y_t = beta1 + (beta2+delta)*x2_t + beta3*x3_t + e_t
# x_2~N(0,1), x_3~N(0,1)
# beta1 = 2, beta2 = 0.5, beta3 = 0.2
# e_t~N(0,1)
# For each sample size, evaluate how power changes with different delta
# simulate the test 500 times 

rm(list=ls(all=T))

# t: sample size, b1: beta1, b2: beta2+delta, b3: beta3
ols_Hypo_type1 <- function(t, delta){
	x2 <- rnorm(t, 0, 1)
	x3 <- rnorm(t, 0, 1)
	b1 <- 2
	b2 <- 0.5 + delta
	b3 <- 0.2
	b  <- t(matrix(c(b1, b2, b3), ncol=3))
	x  <- matrix(c(rep(1,t), x2, x3), ncol=3)
	y  <- x%*%b + rnorm(t, 0, 1)

	coef.lm  <- coef(summary(lm(y~x2+x3)))
	b_hat <- matrix(coef.lm[,1], ncol=1)
	y_hat <- x%*%b_hat
	a <-solve(t(x)%*%x)
	err <-sqrt(t(y-y_hat)%*%(y-y_hat)/(t-3)*diag(a)[2])  # refer to lecture note ch.3: page 56
	t_stat <- (coef.lm[2,1]-0.5)/err
	t_stat
}


# n: simulation times
exe <- function(t, n, delta){
	critical_value <- qt(0.975, t-3)
	test <-numeric(n)
	for (i in 1:n){
		test[i] <- ols_Hypo_type1(t, delta)
	}
	rej <- length(test[abs(test)>critical_value])
	print(rej/n)
}

T <- c(10, 500)    # 2 kinds of sample size
Del <-c(0, 0.1, 0.3, 0.5, 0.7, 1, 1.2) # 7 kinds of delta

y_axis <-seq(from=0, to=1, by=1/(length(Del)-1))
power  <- numeric(length(Del))
plot(Del, y_axis, xlab="Delta", ylab="Power", main="Power Curves (Sample size=10, 500)", type='n')

for (j in 1:length(T)){
	print(j)
	for (i in 1:length(Del)){
		power[i]=exe(T[j], 500, Del[i])
	}
	lines(Del, power, type='b', col=colors()[120+5*(j-1)])
}

text(0.5,0.5, "Blue Curve= sample size 500\n Purple Curve =sample size 10\n", cex=0.8)