Q:

How to deseasonalize quarterly data in R?

I would like to deseasonalize quarterly data with R. I found deseasonalize package, but it says it can deseasonalize only monthly data.
Here's a reproducible example of what I have:
## How to deseasonalize quarterly data? ##
df <- as.data.frame(rnorm(40, 0, 1))
library(lubridate); library(zoo)
row.names(df) <- as.yearqtr(seq(ymd('2000-01-01'), by = '1 quarter', length.out=(40)))
df
seriestobedeseasonalized <- ts(df[,1])

A:

I don't know which filter dimitris_ps is used, but the smoothing degree in his solution is not good. Eventually I devised the following solution with very nice smoothing degree. I used "centered moving averages" as filter. The following code works for any series whose total number of observations is a multiple of 4.
#######deseasonalizeQ: deseasonalize quarterly data #######
# Inspired by excellent work of Jason Delaney on Quarterly Deseasonalize: https://www.youtube.com/watch?v=Jr_2nj6M7L8
# Mutatis mutandis replica of Jason's logic in R 

sales <- ts(c(6,15,10,4,10,18,15,7,14,26,23,12,19,28,25,18,22,34,28,21,24,36,30,20,28,40,35,27))

deseasonalizeQ <- function (x){
x <- ts(x)
#Step1: Centered moving averages: create cma time series having the same length with the original time series x
# cma has 2 NAs on both ends.
cma <- filter(x, filter = c(1/8, 1/4, 1/4, 1/4, 1/8), sides=2)

#Step2: Ratios = Original time series / centered moving averages
ratio <- x/cma

#Step3: Unadjusted 4 seasonal indexes
unadj4si <- ts(1:4)
# floor((length(x)-4)/4)  #"-4" is 4 NA at both ends; below "-1" is due to starting "0:" in multiplication

unadj4si[1] <- mean(ratio[3+4*(0:(floor((length(x)-4)/4) - 1))])
unadj4si[2] <- mean(ratio[4+4*(0:(floor((length(x)-4)/4) - 1))])
unadj4si[3] <- mean(ratio[5+4*(0:(floor((length(x)-4)/4) - 1))])
unadj4si[4] <- mean(ratio[6+4*(0:(floor((length(x)-4)/4) - 1))])

#Step4: Adjusted 4 seasonal indexes
adj4si <- ts(1:4)
adj4si[1] <- unadj4si[1]/mean(c(unadj4si[1],unadj4si[2],unadj4si[3],unadj4si[4]))
adj4si[2] <- unadj4si[2]/mean(c(unadj4si[1],unadj4si[2],unadj4si[3],unadj4si[4]))
adj4si[3] <- unadj4si[3]/mean(c(unadj4si[1],unadj4si[2],unadj4si[3],unadj4si[4]))
adj4si[4] <- unadj4si[4]/mean(c(unadj4si[1],unadj4si[2],unadj4si[3],unadj4si[4]))

#Step5: Propogated adjusted seasonal indexes
propadjsi <- ts(1:length(x))

propadjsi[3+4*(0:(floor((length(x)-4)/4) - 1))] <- adj4si[1]
propadjsi[4+4*(0:(floor((length(x)-4)/4) - 1))] <- adj4si[2]
propadjsi[5+4*(0:(floor((length(x)-4)/4) - 1))] <- adj4si[3]
propadjsi[6+4*(0:(floor((length(x)-4)/4) - 1))] <- adj4si[4]

propadjsi[1] <- adj4si[3]
propadjsi[2] <- adj4si[4]
propadjsi[length(x)-1] <- adj4si[1]
propadjsi[length(x)] <- adj4si[2]

#Step6: Deseasonalized values
out <- x/propadjsi  # deseasonalized = x/propadjsi
out
}

deseasonalizeQ(sales)

#Time Series:Start = 1, End = 28, Frequency = 1  
[1]  6.673117 11.015814  8.941810  6.442787 11.121862 13.218976 13.412714 
[8] 11.274878 15.570607 19.094077 20.566162 19.328362 21.131538 20.562852
[15] 22.354524 28.992543 24.468097 24.969177 25.037067 33.824633 26.692469
[22] 26.437953 26.825429 32.213936 31.141214 29.375503 31.296333 43.488814

###### Plots ########
salesSA <- deseasonalizeQ(sales)
salesSAsalesORJ <- cbind(salesSA, sales)

plot(salesSAsalesORJ, plot.type="single", main="Compare", ylab="values", col=c("blue", "red"), lty=1:2)
legend(10, 40, legend=c("salesSA","sales"), col=c("blue", "red"), lty=1:2)

I checked this deseasonalizeQ with sales(28obs) and growthrate(56obs) time series. Both are OK. Anyway, this solution is long, but its smoothing is considerably better than the above answer.

