#1.) 
library(Rfit)

Y=telephone$calls
x=telephone$year

plot(Y~x,xlab='Rok',ylab='Pocty hovoru')


#(a)
fit<-rfit(Y~x)
summary(fit)

#(b)
summary(lm(Y~x))

plot(Y~x,xlab='Rok',ylab='Pocty hovoru')
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(-1,1))$root
beta_hat


#(d)
Dn=function(b){
R=rank(Y-x*b)
Dn=sum((Y-(x-mean(x))*b)*(R/(n+1)-1/2))
}

beta_hat2=optimize(Dn,c(-1,1))$min
beta_hat2

#(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)

Tn=function(b){
R=rank(Y-x1*b[1]-x2*b[2])
Sn1=sum((x1-mean(x1))*(R/(n+1)-1/2))
Sn2=sum((x2-mean(x2))*(R/(n+1)-1/2))
Tn=sqrt(Sn1^2+Sn2^2)
#Tn=abs(Sn1)+abs(Sn2)
#Tn=max(abs(Sn1),abs(Sn2))
}

optim(c(0,0),Tn)

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

optim(c(0,0),Dn)


#rrs.test(x1,x2,Y)

#(d) H0: beta1=0

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

X=cbind(rep(1,n),x2)
Z_hat=X%*%solve(t(X)%*%X)%*%t(X)%*%x1
Dn_hat=sum((x1-Z_hat)^2)/(n-1)
Tn=sum((x1-Z_hat)*b)/sqrt(n-1)/sqrt(Dn_hat)*sqrt(12)
2*(1-pnorm(Tn))


# aligned test H0: beta1=0

Dn=function(b){
R=rank(Y-(x2-mean(x2))*b)
Dn=sum((Y-(x2-mean(x2))*b)*(R/(n+1)-1/2))
}

beta2_hat=optimize(Dn,c(0,10))$min
e=Y-x2*beta2_hat
R_hat=rank(e)


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


#(e) H0: beta2=0

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

X=cbind(rep(1,n),x1)
Z_hat=X%*%solve(t(X)%*%X)%*%t(X)%*%x2
Dn_hat=sum((x2-Z_hat)^2)/(n-1)
Tn=sum((x2-Dn_hat)*b)/sqrt(n-1)/sqrt(Dn_hat)*sqrt(12)
2*(1-pnorm(abs(Tn)))

#aligned test H0: beta2=0

Dn=function(b){
R=rank(Y-(x1-mean(x1))*b)
Dn=sum((Y-(x1-mean(x1))*b)*(R/(n+1)-1/2))
}

beta1_hat=optimize(Dn,c(0,10))$min
e=Y-x1*beta1_hat
R_hat=rank(e)


Dn=var(x2)
Tn=sum((x2-mean(x2))*R_hat/(n+1))/sqrt(n-1)/sqrt(Dn)*sqrt(12)
2*(1-pnorm(abs(Tn)))