#1.) 
library(Rfit)

#(a)
fit<-rfit(calls~year,data=telephone)
summary(fit)

#(b)
Y=telephone$calls
x=telephone$year

lm(Y~x)

plot(Y~x)
abline(lm(Y~x)$coef)
abline(fit$coef,col=2)

#(c)
n=length(x)
Sn=function(b){
R=rank(Y-x*b)
Sn=sum((x-mean(x))*(R/(n+1)-1/2))
return(Sn)
}

beta_hat=uniroot(Sn,c(0,1))$root
beta_hat


#(d)
J=function(b){
R=rank(Y-(x-mean(x))*b)
J=sum((Y-(x-mean(x))*b)*R)
}

beta_hat2=optimize(J,c(0,10))$min
beta_hat2

Jn=function(b){
R=rank(Y-x*b)
Jn=sum((Y-x*b)*(R/(n+1)-1/2))
return(Jn)
}

beta_hat3=optimize(Jn,c(0,10))$min
beta_hat3

#(e)
e=Y-beta_hat*x
signedrank(e)


########################
#(2)
library(quantreg)

#(a)
Y=bbsalaries$logSalary
x1=bbsalaries$aveWins
x2=bbsalaries$aveGames
n=length(Y)

fit<-rfit(Y~x1+x2)
summary(fit)

#(b)

Ln=function(b){
R=rank(Y-x1*b[1]-x2*b[2])
P1=sum((x1-mean(x1))*R)
P2=sum((x2-mean(x2))*R)
Ln=sqrt(P1^2+P2^2)
#Ln=abs(P1)+abs(P2)
#Ln=max(abs(P1),abs(P2))
}

optim(c(0,0),Ln)

#(c)
J=function(b){
R=rank(Y-x1*b[1]-x2*b[2])
J=sum((Y-(x1-mean(x1))*b[1]-(x2-mean(x2))*b[2])*R)
}

optim(c(0,0),J)


#rrs.test(x1,x2,Y)

# H0: beta2=0

a=rq(Y~x1, tau=-1)
R1=ranks(a)$ranks

X=cbind(rep(1,n),x1)
X2hat=X%*%solve(t(X)%*%X)%*%t(X)%*%x2
Q2=sum((x2-X2hat)^2)/n
T1=sum((x2-X2hat)*R1)/sqrt(n)/sqrt(Q2)*sqrt(12)
2-2*pnorm(abs(T1))

# H0: beta1=0

a=rq(Y~x2, tau=-1)
R2=ranks(a)$ranks

X=cbind(rep(1,n),x2)
X1hat=X%*%solve(t(X)%*%X)%*%t(X)%*%x1
Q2=sum((x1-X1hat)^2)/n
T2=sum((x1-X1hat)*R2)/sqrt(n)/sqrt(Q2)*sqrt(12)
2-2*pnorm(T2)


#aligned testy

# H0: beta2=0

J=function(b){
R=rank(Y-(x1-mean(x1))*b)
J=sum((Y-(x1-mean(x1))*b)*R)
}

beta1=optimize(J,c(0,10))$min
e=Y-x1*beta1
R1=rank(e)


Q=var(x2)
T1=sum((x2-mean(x2))*R1/(n+1))/sqrt(n)/sqrt(Q)*sqrt(12)
2-2*pnorm(abs(T1))


# H0: beta1=0

J=function(b){
R=rank(Y-(x2-mean(x2))*b)
J=sum((Y-(x2-mean(x2))*b)*R)
}

beta2=optimize(J,c(0,10))$min
e=Y-x2*beta2
R2=rank(e)


Q1=var(x1)
T1=sum((x1-mean(x1))*R2/(n+1))/sqrt(n)/sqrt(Q1)*sqrt(12)
2-2*pnorm(abs(T1))





##############################################